Outras formas de gerar um vetor

O operador : é usado para gerar um vetor com todos números que estão entre os operandos e são formados somando números inteiros ao primeiro operando.

1L:10L

##  [1]  1  2  3  4  5  6  7  8  9 10

1.5:9.1

## [1] 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5

A função seq() é usada para criar um vetor de várias formas.

Numa das formas especifica-se o valor inicial, o valor final e o incremento entre elementos do vetor.

seq(1, 9.99, 0.1)

##  [1] 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6
## [18] 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3
## [35] 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0
## [52] 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7
## [69] 7.8 7.9 8.0 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9.0 9.1 9.2 9.3 9.4
## [86] 9.5 9.6 9.7 9.8 9.9

Parâmetros nomeados

Note que chamamos a função passando os parâmetros sem especificação de quais são eles. Eles são recebidos pela função na ordem padrão definida pela função.

Mas no R também é possível passar parâmetros nomeados.

Clique em F1 enquanto tem o cursor em cima da função e veja a ordem dos parâmetros. Veja que outros parâmetros que não utilizamos. Podemos usar length.out ao invés de by:

seq(1L, 10L, length.out = 10L)

##  [1]  1  2  3  4  5  6  7  8  9 10

seq(1L, 10L, length.out = 5L)

## [1]  1.00  3.25  5.50  7.75 10.00

Outro parâmetro, along.with, deixa que criemos um vetor num intervalo determinado e o mesmo número de elementos do vetor passado por este parâmetro.

seq(20, 100, along.with = 1:10)

##  [1]  20.00000  28.88889  37.77778  46.66667  55.55556  64.44444  73.33333
##  [8]  82.22222  91.11111 100.00000

Valores faltantes NA

Valores faltantes ou desconecidos são representados por NA

a <- c(1L,NA)
a

## [1]  1 NA

O valor NA quase sempre contamina os cálculos

media <- mean(a)
media

## [1] NA

mas…

media <- mean(a, na.rm = TRUE)
media

## [1] 1

A exceção são expressões que dão sempre o mesmo resultado independentemente do valor da variável

NA ^ 0

## [1] 1

NA | TRUE

## [1] TRUE

NA & FALSE

## [1] FALSE

A melhor forma de testar se existe um valor NA é is.na

v <- c(1, NA, 2)

is.na(v)

## [1] FALSE  TRUE FALSE

Programação com vetores

As operações do R são vetoriais. Numa operação entre um vetor e um escalar, a operação com o escalar é aplicada a cada elemento do vetor

1:5 * 2

## [1]  2  4  6  8 10

1:10 / 10

##  [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Numa operação com vetores do mesmo tamanho, os elementos são pareados

1:10 * 1:10

##  [1]   1   4   9  16  25  36  49  64  81 100

Programação com vetores - recycling

Outro conceito importante é o de recycling.

Numa operação entre dois vetores de tamanhos diferentes, o vetor menor é repetido ciclicamente de forma a ficar com o mesmo tamanho do vetor maior.

Lembra que toda variável no R é um vetor?

Então… o escalar mostrado no primeiro código do slide anterior é um vetor de 1 elemento que sofre recycling

1:10 * 1:2

##  [1]  1  4  3  8  5 12  7 16  9 20

Estruturas construídas a partir de vetores e listas

Existem estruturas mais complexas na linguagem construídas a partir de vetores e listas.

Data Frame
Matrix
Array
Factor
Estruturas que representam datas
Objetos (no paradigma de orientação a objetos)

Data Frame é a estrutura que mais vamos usar para nossas análises de dados.

Data Frames

Data Frame, e seu primo Tibble, são estruturas muito usadas em análises de dados feitas em R.

O Data Frame consiste em um conjunto de vetores nomeados, com o mesmo número de elementos, que formam uma estrutura retangular, onde cada coluna é um vetor e cada linha n contém o n-ésimo elemento dos vetores.

É similar, em muitas características, a uma tabela de banco de dados.

Essa estrutura é chave no paradigma “Tidy” que usaremos com as bibliotecas Tidyverse

Tibble é uma adaptação do Data Frame para análise de dados. Discutir essas diferenças está fora do escopo do curso. Algumas diferenças serão citadas o longo do material e justificam o uso do Tibble.

df <- 
    data.frame(
        nome = c("João", "Maria", "Zezinho", "Juquinha"), 
        idade = c(7, 8, 9, 10), 
        altura = c(10, 11)
    )
df

##       nome idade altura
## 1     João     7     10
## 2    Maria     8     11
## 3  Zezinho     9     10
## 4 Juquinha    10     11

#tibble não aceita recycling em vetores de tamanho diferente de 1
tib <- 
    #try evita que o erro paralise toda a execução do script
    try(
        tibble(
            nome = c("João", "Maria", "Zezinho", "Juquinha"), 
            idade = c(7, 8, 9, 10), 
            altura = c(10, 11)
        )
    )

## Error : Tibble columns must have consistent lengths, only values of length one are recycled:
## * Length 2: Column `altura`
## * Length 4: Columns `nome`, `idade`
## Backtrace:
##      x
##   1. +-rmarkdown::render("C:/temp/novo/CursoR/Conteudo.Rmd", encoding = "UTF-8")
##   2. | \-knitr::knit(...)
##   3. |   \-knitr:::process_file(text, output)
##   4. |     +-base::withCallingHandlers(...)
##   5. |     +-knitr:::process_group(group)
##   6. |     \-knitr:::process_group.block(group)
##   7. |       \-knitr:::call_block(x)
##   8. |         \-knitr:::block_exec(params)
##   9. |           +-knitr:::in_dir(...)
##  10. |           \-knitr:::evaluate(...)
##  11. |             \-evaluate::evaluate(...)
##  12. |               \-evaluate:::evaluate_call(...)
##  13. |                 +-evaluate:::timing_fn(...)
##  14. |                 +-base:::handle(...)
##  15. |                 +-base::withCallingHandlers(...)
##  16. |                 +-base::withVisible(eval(expr, envir, enclos))
##  17. |                 \-base::eval(expr, envir, enclos)
##  18. |                   \-base::eval(expr, envir, enclos)
##  19. +-base::try(...)
##  20. | \-base::tryCatch(...)
##  21. |   \-base:::tryCatchList(expr, classes, parentenv, handlers)
##  22. |     \-base:::tryCatchOne(expr, names, parentenv, handlers[[1L]])
##  23. |       \-base:::doTryCatch(return(expr), name, parentenv, handler)
##  24. \-tibble::tibble(...)
##  25.   \-tibble:::lst_to_tibble(xlq$output, .rows, .name_repair, lengths = xlq$lengths)
##  26.     \-tibble:::recycle_columns(x, .rows, lengths)

Criando um `tibble` com `tribble()`

É possível criar um Tibble a partir de um código que parece uma tabela, ou seja, criar por linhas e não por colunas.

Isso é feito com a função tribble()

tribble(
  ~nome,       ~idade,      ~altura,
  "Bruno",     41,         1.75,
  "João",      23,         1.80,
  "Maria",     29,         1.70,
  "Zezinho",   31,         1.78
)

## # A tibble: 4 x 3
##   nome    idade altura
##   <chr>   <dbl>  <dbl>
## 1 Bruno      41   1.75
## 2 João       23   1.8 
## 3 Maria      29   1.7 
## 4 Zezinho    31   1.78

Controle de fluxo

A linguagem oferece comandos de controle de fluxo similares aos de outras linguagens.

Podemos dividir os comandos de controle de fluxo em dois tipos:

choices: execução alternativa de comandos
loops: execução repetida de comandos

Choices: `if`, `ifelse`

O comando if funciona para um valor lógico escalar

if (2 + 2 == 4) {
    "2 mais 2 são 4"
} else {
    "2 mais 2 não são 4"
}

## [1] "2 mais 2 são 4"

Note o operador de comparação == e não =

A função if_else (da biblioteca dplyr) funciona de vetorial. if_else é mais rápida que a função ifelse da biblioteca base, mas só aceita argumentos de mesmo tipo no segundo e terceiro parâmetros

jogo_do_pim_silvio_santos <- if_else(
    condition = 1:40 %% 4 == 0 ,
    true =  "PIM",
    false =  as.character(1:40)
)
jogo_do_pim_silvio_santos

##  [1] "1"   "2"   "3"   "PIM" "5"   "6"   "7"   "PIM" "9"   "10"  "11" 
## [12] "PIM" "13"  "14"  "15"  "PIM" "17"  "18"  "19"  "PIM" "21"  "22" 
## [23] "23"  "PIM" "25"  "26"  "27"  "PIM" "29"  "30"  "31"  "PIM" "33" 
## [34] "34"  "35"  "PIM" "37"  "38"  "39"  "PIM"

Note o operador %% e a função de coerção de tipo as.character

Choices: `switch` e `case_when`

A cláusula switch e a função dplyr::case_when evitam que o programador tenha que criar muitos if else aninhados

letra <- "b"

switch(
    letra,
    "a" = "começa com a",
    "b" = "começa com b",
    stop("deu ruim")
)

## [1] "começa com b"

Note que a condição vai sendo testada na ordem e stop gera um erro

case_when serve ao caso vetorial

case_when(
    1:40 %% 10 == 0 ~ "dezena",
    1:40 %% 2 == 0 ~ "par",
    TRUE ~ as.character(1:40)
)

##  [1] "1"      "par"    "3"      "par"    "5"      "par"    "7"     
##  [8] "par"    "9"      "dezena" "11"     "par"    "13"     "par"   
## [15] "15"     "par"    "17"     "par"    "19"     "dezena" "21"    
## [22] "par"    "23"     "par"    "25"     "par"    "27"     "par"   
## [29] "29"     "dezena" "31"     "par"    "33"     "par"    "35"    
## [36] "par"    "37"     "par"    "39"     "dezena"

Loops

A cláusula de loop mais usada e mais versátil é for

for(i in 1:5){ 
    print(i^2)
}

## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25

As cláusulas next e break modificam o comportamento, respectivamente caminhando direto para a próxima iteração e saindo do for

#next vai pra próxima iteração
for(i in 1:5){
    if (i %% 2 == 0){
        next
    }
    print(i)
}

## [1] 1
## [1] 3
## [1] 5

#next sai do loop
for(i in 1:5){
    if (i %% 2 == 0){
        break
    }
    print(i)
}

## [1] 1

Loops: coisa do passado

Vamos ver que quase sempre é desnecessário usar loop para as tarefas que vamos executar.

O caráter vetorial da linguagem, aliado a funcionalidades das bibliotecas, faz com que a grande maioria dos loops sejam desnecessários.

O código fica mais limpo e expressivo e mais rápido. Às vezes MUITO mais rápido. Isso ocorre por motivos além do escopo do curso (alocação de memória, código interpretado x código compilado em C++ etc.)

O código abaixo usa loop e programação funcional, respectivamente. Programação funcional será abordada posteriormente no material.

com_loop <- function(n){
    x <- integer()
    for (i in 1:n){
        x <- c(x, i^2)
    }
    x
}

#programação funcional: aprenderemos posteriomente
sem_loop <- function(n){
    x <- 1:n %>% 
        map_dbl(function(x){x^2})
    x
}

Abaixo as três formas de fazer a mesma conta que terão a performance avaliada

com_loop(5)

## [1]  1  4  9 16 25

sem_loop(5)

## [1]  1  4  9 16 25

(1:5)^2

## [1]  1  4  9 16 25

Loops: coisa do passado (cont.)

A biblioteca bench oferece funções ótimas para avaliar a performance de pedaços pequenos de código.

resultados_perf <- mark(
    sem_loop(1e4),
    com_loop(1e4),
    (1:1e4)^2
)

#aprenderemos o que é %>% e select() posteriormente 
resultados_perf %>% 
    select(expression, min, median, `itr/sec` )

## # A tibble: 3 x 4
##   expression           min   median `itr/sec`
##   <bch:expr>      <bch:tm> <bch:tm>     <dbl>
## 1 sem_loop(10000)   6.65ms   7.56ms    124.  
## 2 com_loop(10000) 118.93ms 121.18ms      5.59
## 3 (1:10000)^2       16.9us   18.5us  25951.

plot(resultados_perf)

Revisão: tipos de vetores

Revisão: operações vetoriais

Qual o valor de v1 ?

v1 <- c(1, 2, 3, 4, 5, 6) * 2

1. (1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6)
1. (2, 4, 6, 8, 10, 12)
1. (2, 2, 3, 4, 5, 6)

Revisão: operações vetoriais

Qual o valor de v2

v2 <-  c(1, 2, 3, 4) + c(0, 1)

1. (1, 3, 3, 5)
1. (1, 3, 3, 4)
1. (1, 2, 3, 4, 0, 1)

Revisão: data frames

Revisão: controle de fluxo

Revisão: if

Complete a lacuna:

if (n_pessoas ____ 0 ){
    print("Não há pessoas")
}
else{
    print("Há pessoas")
}

Revisão: if_else

Complete a lacuna:

x <- 1:10

if_else(
    _________,
    "par",
    "impar"
)

Revisão: if_else

Qual o valor de x ?

times <- c("Flamengo", "Fluminense", "Bahia", "Vasco")

x <- case_when(
    times == "Flamengo" ~ "Flamengo",
    times == "Bahia" ~ "Bahia",
    TRUE ~ "Flumifogo da Gama"
)

Revisão: for

Complete a lacuna para imprimir os quadrados dos números de 1 a 10

for ____{
    print(x^2)
}

Exemplo de simulação: Monty Hall

Monty Hall era uma espécie de Sílvio Santos juvenil (sub 80) americano.

Um dos seus jogos consistia em mostrar três portas ao otár… (ops) convidado. Atrás de uma delas sempre tem um carro. Atrás das outras, alguma brincadeira

Após a escolha inicial, o apresentador revela uma das portas e pergunta se o convidado quer trocar a porta escolhida originalmente.

O que vocês acham? Melhor trocar, manter a escolha original ou tanto faz?

Simulando o Monty Hall

Note o que há de interessante no código (comentado):

A função sample()
A seleção negativa de elementos com -
A função replicate()

set.seed(88)

joga_monty_hall <- function(troca){
    portas <- 1:3
    #sample() sorteia elementos com ou sem reposição
    porta_carro <- sample(portas, size = 1, replace = FALSE)
    primeira_escolha <- 1
    #Seleção negativa (retirando elementos)
    portas_pra_revelar <- portas[-c(porta_carro, primeira_escolha)]
    porta_revelada <- sample( c(portas_pra_revelar, portas_pra_revelar  ), 1)

    if(troca){
        escolha <- portas[-c(primeira_escolha, porta_revelada)]
    }
    else{
        escolha <- primeira_escolha
    }
    
    escolha == porta_carro
        
}

n <- 1000
#replicate executa múltiplas vezes um comando e armazena os resultados em uma estrutura única
troca <- replicate(n = n, joga_monty_hall(troca = TRUE))
fica  <- replicate(n = n, joga_monty_hall(troca = FALSE))

Resultados:

sum(troca)/n

## [1] 0.675

sum(fica)/n

## [1] 0.323

Outra simulação: dá pra passar no CFA sem saber nada?

Vamos ver… mas dá pra simular sem saber quase nada.

Vamos usar uma das funções da família r<familia de distribuição de prob>(). Neste caso, a rbinom, que simula a distribução binomial (aquela que equivale ao evento de jogar n moedas (ou alguma coisa com dois lados) para cima e ver quantas deram cara).

n_simul <- 10000
n_questoes <- 240
min_aprovacao <-  0.6
n_aprovado <- 240 * min_aprovacao
prob_questao <- 0.2

acertos <- rbinom(n = n_simul, size = n_questoes, prob = prob_questao   )

sum(acertos >= n_aprovado)/n_simul

## [1] 0

A chance é praticamente nula.

Na verdade, a grande massa da distribuição fica muito distante.

dado <- enframe(acertos/n_questoes)

mostra_chances <- function(acertos, n_questoes){
    ggplot(enframe(acertos/n_questoes)) +
        geom_density( aes(x = value)) +
        scale_x_continuous(
            labels = percent_format(accuracy = 1), 
            limits = c(0,1),
            breaks = seq(0, 1, 0.1) 
            ) +
        labs(x ="% Acertos") +
        geom_vline(xintercept = min_aprovacao, color = "red") +
        theme_light()
}

mostra_chances(acertos, n_questoes)

Outra simulação: dá pra passar no CFA sabendo a um grau x ?

O exemplo anterior era muito simplista: ninguém chuta tudo.

Imagine que sabemos qual a chance de aparecer uma pergunta onde podemos descartar 0 alternativas, a chance de uma onde descartamos 1 e assim por diante.

Na função rbinom(), a cada simulação, sorteamos uma de duas categorias para cada um dos n eventos. A rmultinom() sorteia uma de k categorias para cada um dos n eventos. Um vetor de k elementos indica a probabilidade de cada categoria.

A rmultinom() nos ajuda a compor provas de acordo com a proficiência de quem vai fazer a prova.

#definindo a chance podermos eliminar 0, 1, 2, ... 4 alternativas
fracao_eliminar_questoes <- c( 0.1, 0.1, 0.2, 0.25 , 0.35 ) 
#definindo o número de questões 
n_questoes_cada_elimina <- t(rmultinom(n_simul, size = n_questoes, fracao_eliminar_questoes))
probs_quando_elimina <- 1/(5:1)
acertos_concatenados <- 
    rbinom( 
        n =  n_simul * 5 , 
        size = as.vector(t(n_questoes_cada_elimina)), 
        prob = probs_quando_elimina  
    )

n_questoes_cada_elimina[1:4,]

##      [,1] [,2] [,3] [,4] [,5]
## [1,]   23   24   50   47   96
## [2,]   24   17   48   64   87
## [3,]   23   30   43   56   88
## [4,]   20   25   51   52   92

acertos_concatenados[1:20]

##  [1]  4  6 17 25 96  4  3 13 34 87  6  5 13 36 88  4  9 17 25 92

Como geramos os acertos concatenados em um vetor, o ideal é transformá-los numa matriz onde cada linha corresponde a uma simulação e cada coluna corresponde aos acertos conseguidos em cada tipo de questão. Usamos a função A matrix(), informando que estamos passando a matriz linha a linha e informando a quantidade de linhas.

matriz_acertos <- matrix(acertos_concatenados, byrow = TRUE, nrow = n_simul )

matriz_acertos[1:5,]

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    4    6   17   25   96
## [2,]    4    3   13   34   87
## [3,]    6    5   13   36   88
## [4,]    4    9   17   25   92
## [5,]    6    5   11   30   92

A rowSums() nos ajuda a ter a quantidade de acertos em cada simulação de prova.

acertos <- rowSums(matriz_acertos)
sum(acertos > n_aprovado)/n_simul

## [1] 0.3131

mostra_chances(acertos, n_questoes)

INTRODUÇÃO À PROGRAMAÇÃO FUNCIONAL COM PURRR

O que é programação funcional?

O fato de que funções são objetos de primeira classe no R, ou seja, objetos que têm propriedades iguais às de qualquer outro, possibilita a progrmação no estilo funcional.

A programação funcional inclui as formas de interação entre vetores e funções a partir de uma função

Neste curso vamos nos concentrar nos Functionals, funções que recebem função como argumento e devolvem um vetor ou uma lista

Função map

A função mais fundamental é map().

Ela recebe uma função e um vetor de \(n\) elementos

A função é chamada para cada elemento do vetor, \(n\) vezes.

Os resultados da aplicação destas \(n\) execuções são devolvidos em uma lista de \(n\) elementos

Exemplo da função map

Uma função é um objeto como outro qualquer e pode ser colocado em uma variável

funcao <- function(x) x^2

map devolve uma lista com o resultado da execução de map em cada elemento

map(.x = 1:10, .f = funcao)

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 4
## 
## [[3]]
## [1] 9
## 
## [[4]]
## [1] 16
## 
## [[5]]
## [1] 25
## 
## [[6]]
## [1] 36
## 
## [[7]]
## [1] 49
## 
## [[8]]
## [1] 64
## 
## [[9]]
## [1] 81
## 
## [[10]]
## [1] 100

`map` pode retornar outros tipos de objeto

map_dbl(.x = 1:10, .f = funcao)

##  [1]   1   4   9  16  25  36  49  64  81 100

map_chr(.x = 1:10, .f = funcao)

##  [1] "1.000000"   "4.000000"   "9.000000"   "16.000000"  "25.000000" 
##  [6] "36.000000"  "49.000000"  "64.000000"  "81.000000"  "100.000000"

map_df(.x = tibble(x = 1:10), .f = funcao)

## # A tibble: 10 x 1
##        x
##    <dbl>
##  1     1
##  2     4
##  3     9
##  4    16
##  5    25
##  6    36
##  7    49
##  8    64
##  9    81
## 10   100

Shortcuts pra funções

Funções podem ser declaradas inline.

map(.x = 1:5, .f = function(x){rnorm(n = 4, mean = x, sd = .2)} )

## [[1]]
## [1] 1.0211304 1.0719202 0.8401777 0.8406449
## 
## [[2]]
## [1] 2.132252 1.703870 1.913016 1.919783
## 
## [[3]]
## [1] 2.731249 2.731211 2.746764 3.003743
## 
## [[4]]
## [1] 3.541141 4.219079 4.103111 3.826659
## 
## [[5]]
## [1] 5.176820 4.772624 4.951106 5.150766

Ou então shortcuts para funções podem ser usados

map(.x = 1:5, .f = ~rnorm(n = 4, mean = .x, sd = .2) )

## [[1]]
## [1] 1.177182 1.243717 1.055257 1.097330
## 
## [[2]]
## [1] 1.979316 1.843530 2.168454 2.225056
## 
## [[3]]
## [1] 2.698754 3.198749 3.099256 3.179662
## 
## [[4]]
## [1] 3.754928 3.706443 4.157722 3.825849
## 
## [[5]]
## [1] 4.799096 4.935156 4.987826 5.313664

Os argumentos são repassados pela map para a função executada. Podemos então passar os argumentos para a map.

map(.x = 1:5, .f = rnorm, n = 4, sd = .2 )

## [[1]]
## [1] 1.163795 1.225138 1.287500 0.721559
## 
## [[2]]
## [1] 2.028696 1.555890 1.915325 2.016978
## 
## [[3]]
## [1] 3.053568 2.672036 2.796469 2.884185
## 
## [[4]]
## [1] 3.867708 3.886027 3.960199 3.927777
## 
## [[5]]
## [1] 4.776481 4.930642 4.631140 4.822234

Versões da `map()` para dois argumentos

A função map()possui versões para mais de um argumento.

map2 e suas modificações para retornos em outros tipos, como map2_dbl(), são preparadas para funções de dois argumentos.

tib_2_arg <- tribble(
  
  ~mean,  ~sd,
  0,       2,
  0,       4,
  1,       2,
  1,       4,
  
)
  

map2(
  .x = tib_2_arg$mean, 
  .y = tib_2_arg$sd, 
  .f = ~rnorm(n = 4, mean = .x, sd = .y) 
)

## [[1]]
## [1]  2.821674 -2.487257 -1.100323 -4.174959
## 
## [[2]]
## [1] -1.1384564 -2.7152693  0.1741787  1.8220994
## 
## [[3]]
## [1] -0.5682192  3.7054536  3.7375878 -1.3283584
## 
## [[4]]
## [1]  3.417435  4.689878  7.298377 -1.680979

Versões da `map()` para mais de dois argumentos

pmap e suas modificações para retornos em outros tipos, como pmap_dbl(), são preparadas um número ilimitado de argumentos.

Uma lista de vetores nomeados deve ser passada para a pmap(). Esses são os argumentos passados para a função.

Lembre que um Data Frame/Tibble é exatamente issouma lista de vetores nomeados.

O detalhe é que os nomes dos vetores deve bater com o nome dos parâmetros da função

tib_n_arg <- tribble(
  
  ~mean,  ~sd,  ~n,
  0,       2,   5,
  0,       4,   5,
  1,       2,   5,
  1,       4,   5,
  0,       2,   10,
  0,       4,   10,
  1,       2,   10,
  1,       4,   10
  
)

pmap(.l = tib_n_arg, .f = rnorm)

## [[1]]
## [1]  0.8674085 -4.5345372 -0.9266187  2.8163247 -1.0953331
## 
## [[2]]
## [1]  5.136385 -3.057098 -2.295600 -1.794702 -1.311037
## 
## [[3]]
## [1]  1.532084 -7.068620  2.369469 -1.790043 -2.132151
## 
## [[4]]
## [1]  5.5805530  0.8476443  1.8301349  3.0455240 -6.1892686
## 
## [[5]]
##  [1]  2.5492313  0.7512555 -2.1446984  3.1894252 -3.5840746  2.2253177
##  [7] -0.6031117  0.6509431  0.6480136  2.1763102
## 
## [[6]]
##  [1]  1.2978069 -1.7907063 -3.2601022  4.0173588 -2.3454343 -1.0336582
##  [7]  0.6735530 -2.6442484  0.5021318 -2.2966865
## 
## [[7]]
##  [1]  3.2827845  0.8395007 -0.7688994  0.3239543  1.9593301  4.4368960
##  [7] -3.6475308  1.1799538  1.9458926  3.2373521
## 
## [[8]]
##  [1]  2.330355  5.263333  1.419372  0.371877 -6.420001  5.127040 -5.708234
##  [8]  2.881843  2.085133  3.274547

Vantagens da programação funcional

A programação funcional deixa o código mais expressivo e maioria das vezes mais performático

Existem algumas características avançadas da programação funcional que facilitam bastante o trabalho.

A principal delas é a facilidade com que podemos converter um código que roda em série para um código que roda em paralelocom a biblioteca furrr, que veremos adiante.

TIDY DATA (obtenção e organização dos dados)

Um habilidade subestimada

Dentro de todo o hype envolvendo Data Science, surgem as buzz words mais mirabolantes: machine learning, AI, Deep Learning…

Tudo isso é legal, mas a habilidade de preparar os dados para os modelos, preparar os hiperparâmetros e especificações alternativas ainda melhora muito a análise. Anote mais uma buzz word: FEATURE ENGINEERING

A agilidade de se tentar abordagens alterativas com os dados cresce muito quando dominamos a arte de manipular os data frames. Por isso o peso grande dado neste curso.

Organizando os dados de forma tidy

Arrumar os dados de forma que as linhas sejam eventos e as colunas sejam atributos do evento ajuda muito a rodar modelos e construir visualizações eficientemente. Esta forma de organização foi chamada de Tidy data por Hadley Wickham, o criador do conjunto de bibliotecas Tidyverse

O que é o evento e o que é o atributo pode variar até para diferentes usos do mesmo dado. Mas a prática ajuda a determinar isso.

Tratamento de dados em passos: operador Pipe (`%>%`)

O operador pipe, representado por %>% é extremamente útil para análise de dados no R com uso das bibliotecas Tidyverse

Normalmente os tratamentos de dados são feitos em múltiplos passos encadeados:

#dados de exemplo
head(gapminder)

## # A tibble: 6 x 6
##   country     continent  year lifeExp      pop gdpPercap
##   <fct>       <fct>     <int>   <dbl>    <int>     <dbl>
## 1 Afghanistan Asia       1952    28.8  8425333      779.
## 2 Afghanistan Asia       1957    30.3  9240934      821.
## 3 Afghanistan Asia       1962    32.0 10267083      853.
## 4 Afghanistan Asia       1967    34.0 11537966      836.
## 5 Afghanistan Asia       1972    36.1 13079460      740.
## 6 Afghanistan Asia       1977    38.4 14880372      786.

Vamos imaginar que queremos a média de PIB per capita por continente em 2007.

Note quanto código desnecessário há nestas linhas: variáveis que não precisavam ser nomeadas nem passadas explicitamente como parâmetro.

Este código desnecessário causa fadiga no programador, confunde o próprio autor do código e confunde ainda mais o leitor posterior do código.

#vamos cobrir essas funções de tratamento posteriormente
gapminder_07 <- filter(gapminder, year == 2007)
gapminder_07_group_continente <- group_by(gapminder_07, continent)
gapminder_media_gdp_continente <- summarise(
    gapminder_07_group_continente, media_gdp = sum(gdpPercap * pop)/sum(pop)
)
resultado <- arrange(gapminder_media_gdp_continente, desc(media_gdp))

resultado

## # A tibble: 5 x 2
##   continent media_gdp
##   <fct>         <dbl>
## 1 Oceania      32885.
## 2 Europe       25244.
## 3 Americas     21603.
## 4 Asia          5432.
## 5 Africa        2561.

Tratamento de dados em passos: operador Pipe (`%>%`) (cont.)

O operador pipe %>% faz o seguinte:

x %>% y(z) = y(x,z)

Ou seja, o primeiro operando é enfiado como primeiro parâmetro da função que está no segundo operando.

Isso faz com que possamos escrever o código anterior assim:

resultado <- gapminder %>% 
    filter(year == 2007) %>% 
    group_by(continent) %>% 
    summarise(
        media_gdp = sum(gdpPercap * pop) / sum(pop)
    ) %>% 
    arrange(desc(media_gdp))
    
resultado

## # A tibble: 5 x 2
##   continent media_gdp
##   <fct>         <dbl>
## 1 Oceania      32885.
## 2 Europe       25244.
## 3 Americas     21603.
## 4 Asia          5432.
## 5 Africa        2561.

Note que agora podemos interpretar o código facilmente como uma série de comandos de tratamento em cima dos dados.

Não é por coincidência que as funções de tratamento das bibliotecas tidyverse que veremos adiante são verbos e recebem os dados como primeiro parâmetro.

Agora o mais importante de tudo: O ATALHO PARA O %>% É CTRL + SHIFT + M

Um mapa conceitual da bibioteca de transformação de dados `dplyr`

dplyr é a biblioteca mais básica do conjunto tidyverse

Fonte: (Courtiol 2019)

CRAN: uma Disneylândia dos dados?

CRAN é o repositório de bibliotecas mantido pelo R com contribuição de populares.

Além de funcionalidades estatísticas e funcionalidades para lidar com dados, há dados e funcionalidades para buscar dados online.

Usaremos várias das bases como exemplo.

A primeira é a do Banco Mundial, muito rica para quem gosta de dados socioeconômicos: wbstats

Para acessar um indicador precisamos achá-lo na base de indicadores com a função wbsearch()

#pattern é uma expressão regular. \\ serve para dizer que "(" é mesmo "(" 
#e não o ( usado nas operações de expressão regular (fora do escopo do curso)
indicadores <- wbsearch(pattern = "GINI index \\(World Bank estimate\\)")

indicadores

##      indicatorID                        indicator
## 1348 SI.POV.GINI GINI index (World Bank estimate)

Sabendo o ID do indicador, podemos consultá-lo com a função wb()

#mrv é most recent values. Pode ser usado para buscar os n valores mais recentes
gini = wb(indicator = "SI.POV.GINI", mrv= 10, POSIXct = TRUE)

head(gini)

##     iso3c date value indicatorID                        indicator iso2c
## 486   ALB 2012  29.0 SI.POV.GINI GINI index (World Bank estimate)    AL
## 490   ALB 2008  30.0 SI.POV.GINI GINI index (World Bank estimate)    AL
## 497   DZA 2011  27.6 SI.POV.GINI GINI index (World Bank estimate)    DZ
## 530   AGO 2008  42.7 SI.POV.GINI GINI index (World Bank estimate)    AO
## 541   ARG 2017  41.2 SI.POV.GINI GINI index (World Bank estimate)    AR
## 542   ARG 2016  42.0 SI.POV.GINI GINI index (World Bank estimate)    AR
##       country    date_ct granularity
## 486   Albania 2012-01-01      annual
## 490   Albania 2008-01-01      annual
## 497   Algeria 2011-01-01      annual
## 530    Angola 2008-01-01      annual
## 541 Argentina 2017-01-01      annual
## 542 Argentina 2016-01-01      annual

dplyr: Modificar -> Colunas -> Nomes e posições

Funções básicas de tratamento (dplyr): select()

A função select() é usada para selecionar colunas do data frame/tibble

glimpse(gini)

## Observations: 682
## Variables: 9
## $ iso3c       <chr> "ALB", "ALB", "DZA", "AGO", "ARG", "ARG", "ARG", "...
## $ date        <chr> "2012", "2008", "2011", "2008", "2017", "2016", "2...
## $ value       <dbl> 29.0, 30.0, 27.6, 42.7, 41.2, 42.0, 41.7, 41.0, 41...
## $ indicatorID <chr> "SI.POV.GINI", "SI.POV.GINI", "SI.POV.GINI", "SI.P...
## $ indicator   <chr> "GINI index (World Bank estimate)", "GINI index (W...
## $ iso2c       <chr> "AL", "AL", "DZ", "AO", "AR", "AR", "AR", "AR", "A...
## $ country     <chr> "Albania", "Albania", "Algeria", "Angola", "Argent...
## $ date_ct     <date> 2012-01-01, 2008-01-01, 2011-01-01, 2008-01-01, 2...
## $ granularity <chr> "annual", "annual", "annual", "annual", "annual", ...

gini_select <- gini %>% 
    select(country, date, value, iso3c)

head(gini_select)

##       country date value iso3c
## 486   Albania 2012  29.0   ALB
## 490   Albania 2008  30.0   ALB
## 497   Algeria 2011  27.6   DZA
## 530    Angola 2008  42.7   AGO
## 541 Argentina 2017  41.2   ARG
## 542 Argentina 2016  42.0   ARG

É possível usar a seleção negativa assim como fizemos com vetores

gini_select2 <- gini_select %>% 
    select(-iso3c)

head(gini_select2)

##       country date value
## 486   Albania 2012  29.0
## 490   Albania 2008  30.0
## 497   Algeria 2011  27.6
## 530    Angola 2008  42.7
## 541 Argentina 2017  41.2
## 542 Argentina 2016  42.0

Funções básicas de tratamento (dplyr): select() (cont.)

Algumas funções helpers da dplyr nos ajudam a usar a função select e são muito úteis para tratamentos mais elaborados.

Pra mostrar mais funcionalidades da função select, vamos usar uma base com dados eleitorais brasileiros, que retorna mais colunas

candidatos <- candidate_fed(2018)

## Processing the data...

## Done.

glimpse(candidatos)

## Observations: 29,113
## Variables: 58
## $ DATA_GERACAO                   <chr> "02/03/2020", "02/03/2020", "02...
## $ HORA_GERACAO                   <time> 19:10:38, 19:10:38, 19:10:38, ...
## $ ANO_ELEICAO                    <dbl> 2018, 2018, 2018, 2018, 2018, 2...
## $ COD_TIPO_ELEICAO               <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ NOME_TIPO_ELEICAO              <chr> "ELEIÇÃO ORDINÁRIA", "ELEIÇÃO O...
## $ NUM_TURNO                      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ COD_ELEICAO                    <dbl> 297, 297, 297, 297, 297, 297, 2...
## $ DESCRICAO_ELEICAO              <chr> "Eleições Gerais Estaduais 2018...
## $ DATA_ELEICAO                   <chr> "07/10/2018", "07/10/2018", "07...
## $ ABRANGENCIA                    <chr> "ESTADUAL", "ESTADUAL", "ESTADU...
## $ SIGLA_UF                       <chr> "AC", "AC", "AC", "AC", "AC", "...
## $ SIGLA_UE                       <chr> "AC", "AC", "AC", "AC", "AC", "...
## $ DESCRICAO_UE                   <chr> "ACRE", "ACRE", "ACRE", "ACRE",...
## $ CODIGO_CARGO                   <dbl> 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7...
## $ DESCRICAO_CARGO                <chr> "DEPUTADO FEDERAL", "DEPUTADO F...
## $ SEQUENCIAL_CANDIDATO           <dbl> 10000600848, 10000600844, 10000...
## $ NUMERO_CANDIDATO               <dbl> 1010, 5012, 10555, 28190, 23555...
## $ NOME_CANDIDATO                 <chr> "MANUEL MARCOS CARVALHO DE MESQ...
## $ NOME_URNA_CANDIDATO            <chr> "PASTOR MANUEL MARCOS", "GLÓRIA...
## $ NOME_SOCIAL_CANDIDATO          <chr> "#NULO#", "#NULO#", "#NULO#", "...
## $ CPF_CANDIDATO                  <chr> "36089427268", "02319154205", "...
## $ EMAIL_CANDIDATO                <chr> "MANUELMARCOS323@GMAIL.COM", "G...
## $ COD_SITUACAO_CANDIDATURA       <dbl> 12, 12, 12, 12, 12, 12, 12, 12,...
## $ DES_SITUACAO_CANDIDATURA       <chr> "APTO", "APTO", "APTO", "APTO",...
## $ COD_DETALHE_SITUACAO_CAND      <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ DES_DETALHE_SITUACAO_CAND      <chr> "DEFERIDO", "DEFERIDO", "DEFERI...
## $ TIPO_AGREMIACAO                <chr> "COLIGAÇÃO", "COLIGAÇÃO", "COLI...
## $ NUMERO_PARTIDO                 <dbl> 10, 50, 10, 28, 23, 10, 11, 77,...
## $ SIGLA_PARTIDO                  <chr> "PRB", "PSOL", "PRB", "PRTB", "...
## $ NOME_PARTIDO                   <chr> "PARTIDO REPUBLICANO BRASILEIRO...
## $ CODIGO_LEGENDA                 <dbl> 10000050031, 10000050031, 10000...
## $ NOME_COLIGACAO                 <chr> "CHAPINHA - FPA II", "CHAPINHA ...
## $ COMPOSICAO_LEGENDA             <chr> "PSOL / PV / PPL / PRP / PRB / ...
## $ CODIGO_NACIONALIDADE           <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ DESCRICAO_NACIONALIDADE        <chr> "BRASILEIRA NATA", "BRASILEIRA ...
## $ SIGLA_UF_NASCIMENTO            <chr> "CE", "AC", "AC", "AC", "AC", "...
## $ CODIGO_MUNICIPIO_NASCIMENTO    <dbl> -3, -3, -3, -3, -3, -3, -3, -3,...
## $ NOME_MUNICIPIO_NASCIMENTO      <chr> "TIANGUA", "RIO BRANCO", "RIO B...
## $ DATA_NASCIMENTO                <chr> "07/05/1972", "16/02/1995", "16...
## $ IDADE_DATA_POSSE               <dbl> 46, 23, 33, 35, 49, 43, 52, 47,...
## $ NUM_TITULO_ELEITORAL_CANDIDATO <chr> "010751102267", "006492162496",...
## $ CODIGO_SEXO                    <dbl> 2, 4, 2, 2, 2, 2, 2, 4, 2, 2, 2...
## $ DESCRICAO_SEXO                 <chr> "MASCULINO", "FEMININO", "MASCU...
## $ COD_GRAU_INSTRUCAO             <dbl> 6, 6, 6, 6, 6, 8, 6, 6, 4, 5, 8...
## $ DESCRICAO_GRAU_INSTRUCAO       <chr> "ENSINO MÉDIO COMPLETO", "ENSIN...
## $ CODIGO_ESTADO_CIVIL            <dbl> 3, 1, 1, 1, 1, 3, 3, 3, 1, 1, 9...
## $ DESCRICAO_ESTADO_CIVIL         <chr> "CASADO(A)", "SOLTEIRO(A)", "SO...
## $ CODIGO_COR_RACA                <chr> "03", "01", "03", "03", "03", "...
## $ DESCRICAO_COR_RACA             <chr> "PARDA", "BRANCA", "PARDA", "PA...
## $ CODIGO_OCUPACAO                <dbl> 278, 999, 541, 601, 601, 277, 3...
## $ DESCRICAO_OCUPACAO             <chr> "VEREADOR", "OUTROS", "MECÂNICO...
## $ DESPESA_MAX_CAMPANHA           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, ...
## $ COD_SIT_TOT_TURNO              <dbl> 3, 5, 5, 4, 5, 5, 5, 5, 5, 4, 5...
## $ DESC_SIT_TOT_TURNO             <chr> "ELEITO POR MÉDIA", "SUPLENTE",...
## $ SITUACAO_REELEICAO             <chr> "N", "N", "N", "N", "N", "S", "...
## $ SITUACAO_DECLARAR_BENS         <chr> "N", "N", "S", "N", "S", "S", "...
## $ NUMERO_PROTOCOLO_CANDIDATURA   <dbl> -1, -1, -1, -1, -1, -1, -1, -1,...
## $ NUMERO_PROCESSO                <chr> "06002849320186010000", "060028...

Funções básicas de tratamento (dplyr): select() - helpers

candidatos_select <- candidatos %>% 
    select(starts_with("NOME"))

datatable(head(candidatos_select))

candidatos_select <- candidatos %>% 
    select(ends_with("candidato"))

datatable(head(candidatos_select))

candidatos_select <- candidatos %>% 
    select(contains("municipio"))

datatable(head(candidatos_select))

candidatos_select <- candidatos %>% 
    select(ends_with("candidato"))

datatable(head(candidatos_select))

Funções básicas de tratamento (dplyr): select() - helpers (cont.)

A função helper num_range ajuda a encontrar colunas do tipo prefixo_n. Isso é muito comum em bases de dados

A biblioteca worldmet retorna dados de estações meteorológicas espalhadas pelo planeta

Primeiro é necessário encontrar o código da base desejada

estacao <- getMeta("heathrow", returnMap = TRUE)

estacao

Funções básicas de tratamento (dplyr): select() - helpers (cont.)

A função abaixo retorna os dados de uma estação. Veja que alguns campos têm um sufixo _n

dados_heathrow <- importNOAA(code = "037720-99999", year = 2019,
precip = TRUE, PWC = FALSE, parallel = TRUE)


glimpse(dados_heathrow)

## Observations: 8,760
## Variables: 26
## $ date        <dttm> 2019-01-01 00:00:00, 2019-01-01 01:00:00, 2019-01...
## $ usaf        <chr> "037720", "037720", "037720", "037720", "037720", ...
## $ wban        <chr> "99999", "99999", "99999", "99999", "99999", "9999...
## $ code        <chr> "037720-99999", "037720-99999", "037720-99999", "0...
## $ station     <chr> "HEATHROW", "HEATHROW", "HEATHROW", "HEATHROW", "H...
## $ lat         <dbl> 51.47967, 51.47967, 51.47967, 51.47967, 51.47967, ...
## $ lon         <dbl> -0.4573333, -0.4573333, -0.4573333, -0.4573333, -0...
## $ elev        <dbl> 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25...
## $ wd          <dbl> 283.706052, 286.994553, 290.000000, 290.000000, 27...
## $ ws          <dbl> 3.266667, 3.433333, 4.100000, 3.433333, 3.266667, ...
## $ ceil_hgt    <dbl> 657.3333, 667.3333, 728.0000, 768.0000, 778.3333, ...
## $ visibility  <dbl> 22000, 28000, 45000, 45000, 35000, 35000, 28000, 2...
## $ air_temp    <dbl> 8.666667, 8.933333, 9.000000, 9.000000, 7.966667, ...
## $ dew_point   <dbl> 4.9333333, 4.9666667, 4.7666667, 4.8000000, 4.8000...
## $ atmos_pres  <dbl> 1034.8, 1034.5, 1034.6, 1034.5, 1034.4, 1033.9, 10...
## $ RH          <dbl> 77.67802, 76.42835, 75.06237, 75.23079, 80.81974, ...
## $ cl_1        <dbl> 7.666667, 8.000000, 7.666667, 8.000000, 7.666667, ...
## $ cl_1_height <dbl> 657.3333, 667.3333, 728.0000, 768.0000, 778.3333, ...
## $ cl_2        <dbl> 8.000000, NA, 8.000000, NA, NA, NA, 7.000000, 7.00...
## $ cl_2_height <dbl> 792, NA, 914, NA, NA, NA, 1006, 1036, NA, 720, 121...
## $ cl_3        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ cl_3_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ precip_12   <dbl> NA, NA, NA, NA, NA, NA, 0, NA, NA, NA, NA, NA, NA,...
## $ precip_6    <dbl> 0, NA, NA, NA, NA, NA, 0, NA, NA, NA, NA, NA, 0, N...
## $ cl          <dbl> 8.000000, 8.000000, 8.000000, 8.000000, 7.666667, ...
## $ precip      <dbl> NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

A função helper num_range ajuda a selecionar essas colunas com prefixo comum e um sufixo numérico

dados_heathrow_select <- dados_heathrow %>% 
    select( 
        date, 
        num_range("cl_", 1:3 ), 
        num_range("precip_", c(6, 12))  
    )


head(dados_heathrow_select)

## # A tibble: 6 x 6
##   date                 cl_1  cl_2  cl_3 precip_6 precip_12
##   <dttm>              <dbl> <dbl> <dbl>    <dbl>     <dbl>
## 1 2019-01-01 00:00:00  7.67     8    NA        0        NA
## 2 2019-01-01 01:00:00  8       NA    NA       NA        NA
## 3 2019-01-01 02:00:00  7.67     8    NA       NA        NA
## 4 2019-01-01 03:00:00  8       NA    NA       NA        NA
## 5 2019-01-01 04:00:00  7.67    NA    NA       NA        NA
## 6 2019-01-01 05:00:00  6       NA    NA       NA        NA

Outra função útil é a everything(), que ajuda, por exemplo, a passar algumas colunas para o início do tibble.

dados_heathrow_select <- dados_heathrow %>% 
    select( 
        date, 
        air_temp,
        everything() 
    )


glimpse(dados_heathrow_select)

## Observations: 8,760
## Variables: 26
## $ date        <dttm> 2019-01-01 00:00:00, 2019-01-01 01:00:00, 2019-01...
## $ air_temp    <dbl> 8.666667, 8.933333, 9.000000, 9.000000, 7.966667, ...
## $ usaf        <chr> "037720", "037720", "037720", "037720", "037720", ...
## $ wban        <chr> "99999", "99999", "99999", "99999", "99999", "9999...
## $ code        <chr> "037720-99999", "037720-99999", "037720-99999", "0...
## $ station     <chr> "HEATHROW", "HEATHROW", "HEATHROW", "HEATHROW", "H...
## $ lat         <dbl> 51.47967, 51.47967, 51.47967, 51.47967, 51.47967, ...
## $ lon         <dbl> -0.4573333, -0.4573333, -0.4573333, -0.4573333, -0...
## $ elev        <dbl> 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25...
## $ wd          <dbl> 283.706052, 286.994553, 290.000000, 290.000000, 27...
## $ ws          <dbl> 3.266667, 3.433333, 4.100000, 3.433333, 3.266667, ...
## $ ceil_hgt    <dbl> 657.3333, 667.3333, 728.0000, 768.0000, 778.3333, ...
## $ visibility  <dbl> 22000, 28000, 45000, 45000, 35000, 35000, 28000, 2...
## $ dew_point   <dbl> 4.9333333, 4.9666667, 4.7666667, 4.8000000, 4.8000...
## $ atmos_pres  <dbl> 1034.8, 1034.5, 1034.6, 1034.5, 1034.4, 1033.9, 10...
## $ RH          <dbl> 77.67802, 76.42835, 75.06237, 75.23079, 80.81974, ...
## $ cl_1        <dbl> 7.666667, 8.000000, 7.666667, 8.000000, 7.666667, ...
## $ cl_1_height <dbl> 657.3333, 667.3333, 728.0000, 768.0000, 778.3333, ...
## $ cl_2        <dbl> 8.000000, NA, 8.000000, NA, NA, NA, 7.000000, 7.00...
## $ cl_2_height <dbl> 792, NA, 914, NA, NA, NA, 1006, 1036, NA, 720, 121...
## $ cl_3        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ cl_3_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ precip_12   <dbl> NA, NA, NA, NA, NA, NA, 0, NA, NA, NA, NA, NA, NA,...
## $ precip_6    <dbl> 0, NA, NA, NA, NA, NA, 0, NA, NA, NA, NA, NA, 0, N...
## $ cl          <dbl> 8.000000, 8.000000, 8.000000, 8.000000, 7.666667, ...
## $ precip      <dbl> NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

Funções básicas de tratamento (dplyr): `rename()`

rename() é usada para modificar os nomes das colunas. Ela renomeia as colunas indicadas e mantném as outras.

Funções básicas de tratamento (dplyr): `rename()`

dados_heathrow %>% 
    rename(
        data = date,
        estacao = station,
        temperatura = air_temp
    ) %>% 
    head()

## # A tibble: 6 x 26
##   data                usaf  wban  code  estacao   lat    lon  elev    wd
##   <dttm>              <chr> <chr> <chr> <chr>   <dbl>  <dbl> <dbl> <dbl>
## 1 2019-01-01 00:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  284.
## 2 2019-01-01 01:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  287.
## 3 2019-01-01 02:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  290 
## 4 2019-01-01 03:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  290 
## 5 2019-01-01 04:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  276.
## 6 2019-01-01 05:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  267.
## # ... with 17 more variables: ws <dbl>, ceil_hgt <dbl>, visibility <dbl>,
## #   temperatura <dbl>, dew_point <dbl>, atmos_pres <dbl>, RH <dbl>,
## #   cl_1 <dbl>, cl_1_height <dbl>, cl_2 <dbl>, cl_2_height <dbl>,
## #   cl_3 <dbl>, cl_3_height <dbl>, precip_12 <dbl>, precip_6 <dbl>,
## #   cl <dbl>, precip <dbl>

Funções básicas de tratamento (dplyr): `rename()` x `select()`

select() também pode ser usado para renomear colunas, mas mantém apenas as colunas citadas

dados_heathrow %>% 
    select(
        data = date,
        estacao = station,
        temperatura = air_temp
    ) %>% 
    head()

## # A tibble: 6 x 3
##   data                estacao  temperatura
##   <dttm>              <chr>          <dbl>
## 1 2019-01-01 00:00:00 HEATHROW        8.67
## 2 2019-01-01 01:00:00 HEATHROW        8.93
## 3 2019-01-01 02:00:00 HEATHROW        9   
## 4 2019-01-01 03:00:00 HEATHROW        9   
## 5 2019-01-01 04:00:00 HEATHROW        7.97
## 6 2019-01-01 05:00:00 HEATHROW        6.6

dplyr: Modificar -> Colunas -> Valores

A função mutate() é usada para criar novas colunas no tibble, mantendo as outras.

Funções básicas de tratamento (dplyr): `mutate()`

Notando que a coluna DATA_ELEICAO é um caracter, vamos criar uma coluna de tipo data.

typeof(candidatos$DATA_ELEICAO)

## [1] "character"

O jeito mais fácil de fazer isso é usando uma das funções da biblioteca lubridate

candidatos_com_data <- candidatos %>% 
    mutate(DATA_ELEICAO_TIPO_DATA = dmy(DATA_ELEICAO)) %>% 
    select(DATA_ELEICAO, DATA_ELEICAO_TIPO_DATA)

head(candidatos_com_data)

## # A tibble: 6 x 2
##   DATA_ELEICAO DATA_ELEICAO_TIPO_DATA
##   <chr>        <date>                
## 1 07/10/2018   2018-10-07            
## 2 07/10/2018   2018-10-07            
## 3 07/10/2018   2018-10-07            
## 4 07/10/2018   2018-10-07            
## 5 07/10/2018   2018-10-07            
## 6 07/10/2018   2018-10-07

É possível substituir um campo existente

candidatos_com_data <- candidatos %>% 
    mutate(DATA_ELEICAO = dmy(DATA_ELEICAO)) %>% 
    select(DATA_ELEICAO)

head(candidatos_com_data)

## # A tibble: 6 x 1
##   DATA_ELEICAO
##   <date>      
## 1 2018-10-07  
## 2 2018-10-07  
## 3 2018-10-07  
## 4 2018-10-07  
## 5 2018-10-07  
## 6 2018-10-07

Funções básicas de tratamento (dplyr): `mutate()`

funções derivadas da mutate possibilitam a alteração de várias colunas ao mesmo tempo, usando os mesmos helpers que já vimos para a select e uma função à escolha

candidatos_com_data <- candidatos %>% 
    mutate_at(vars(starts_with("DATA_")), dmy ) %>% 
    select(starts_with("DATA_"))


head(candidatos_com_data)

## # A tibble: 6 x 3
##   DATA_GERACAO DATA_ELEICAO DATA_NASCIMENTO
##   <date>       <date>       <date>         
## 1 2020-03-02   2018-10-07   1972-05-07     
## 2 2020-03-02   2018-10-07   1995-02-16     
## 3 2020-03-02   2018-10-07   1985-06-16     
## 4 2020-03-02   2018-10-07   1983-02-19     
## 5 2020-03-02   2018-10-07   1969-07-28     
## 6 2020-03-02   2018-10-07   1975-10-20

Funções básicas de tratamento (dplyr) `mutate()` (cont.):

Outras funções úteis são as que fazem operações acumuladas e as operações de lag() e lead().

BETS é uma biblioteca criada pela FGV que dá acesso a séries temporais econômicas

series <- BETS::BETSsearch("exchange dollar")

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

## 
## BETS-package: Found 33 out of 18706 time series.

series

## # A tibble: 33 x 7
##    code  description           unit   periodicity start  last_value source 
##    <chr> <chr>                 <chr>  <chr>       <chr>  <chr>      <chr>  
##  1 1     Exchange rate - Free~ c.m.u~ D           28/11~ 26/03/2018 Sisbac~
##  2 10813 Exchange rate - Free~ c.m.u~ D           28/11~ 03/05/2018 Sisbac~
##  3 11753 Real effective excha~ Index  M           31/01~ mar/2018   BCB-DS~
##  4 11758 Real effective excha~ Index  M           31/01~ mar/2018   BCB-DS~
##  5 11763 Real effective excha~ Index  M           31/01~ mar/2018   BCB-DS~
##  6 11768 Real effective excha~ Index  M           31/01~ mar/2018   BCB-DS~
##  7 17887 Exchange rate (perio~ US$/c~ M           31/01~ nov/2016   IMF-IFS
##  8 17888 Exchange rate (perio~ US$/c~ M           31/01~ nov/2016   IMF-IFS
##  9 18441 Exchange rate (end o~ US$/c~ M           31/01~ nov/2016   IMF-IFS
## 10 18442 Exchange rate (end o~ US$/c~ M           31/01~ nov/2016   IMF-IFS
## # ... with 23 more rows

No código abaixo, calculamos o retorno da série, a volatilidade histórica e a volatilidade EWMA

dolar <- BETS::BETSget(1) 

dolar_com_vol <- dolar %>% 
    filter(date > ymd("1994-07-01")) %>% 
    arrange(date) %>% 
    mutate(
        retorno = (value - lag(value))/value,
        retorno_quad = retorno^2,
        dia = row_number(),
        fator_ewma = (1/0.94)^dia*1e-20
    ) %>% 
    filter(!is.na(retorno)) %>% 
    mutate(vol = sqrt(cumvar(retorno)) * sqrt(252) ) %>% 
    mutate(vol_ewma = sqrt(cumsum(retorno_quad * fator_ewma)/cumsum(fator_ewma)) * sqrt(252) ) %>% 
    rename(dolar = value) %>% 
    select(
        date,
        dolar,
        retorno,
        vol,
        vol_ewma
    )



datatable(dolar_com_vol) %>% 
    formatPercentage(c("retorno", "vol", "vol_ewma"), 2)

dolar_ajeitado <- dolar_com_vol %>% 
    gather(variavel, valor, - date)


dolar_ajeitado %>% 
    ggplot() +
    geom_line(aes(x = date, y = valor)) +
    facet_grid( variavel ~ . , scales = "free") +
    theme_light()

Funções básicas de tratamento (dplyr) `transmute()`:

transmute() cria colunas e mantém apenas as colunas citadas

gapminder %>% 
  transmute(
    ano = year,
    pais = country,
    pib = gdpPercap * pop
  ) %>% 
  head()

## # A tibble: 6 x 3
##     ano pais                 pib
##   <int> <fct>              <dbl>
## 1  1952 Afghanistan  6567086330.
## 2  1957 Afghanistan  7585448670.
## 3  1962 Afghanistan  8758855797.
## 4  1967 Afghanistan  9648014150.
## 5  1972 Afghanistan  9678553274.
## 6  1977 Afghanistan 11697659231.

dplyr: Modificar -> Linhas -> Posição

A função arrange serve para ordenar as linhas do tibble.

Funções básicas de tratamento (dplyr) `arrange()`:

dados_ordenados <- dados_heathrow %>% 
    arrange(date)

head(dados_ordenados)

## # A tibble: 6 x 26
##   date                usaf  wban  code  station   lat    lon  elev    wd
##   <dttm>              <chr> <chr> <chr> <chr>   <dbl>  <dbl> <dbl> <dbl>
## 1 2019-01-01 00:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  284.
## 2 2019-01-01 01:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  287.
## 3 2019-01-01 02:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  290 
## 4 2019-01-01 03:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  290 
## 5 2019-01-01 04:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  276.
## 6 2019-01-01 05:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  267.
## # ... with 17 more variables: ws <dbl>, ceil_hgt <dbl>, visibility <dbl>,
## #   air_temp <dbl>, dew_point <dbl>, atmos_pres <dbl>, RH <dbl>,
## #   cl_1 <dbl>, cl_1_height <dbl>, cl_2 <dbl>, cl_2_height <dbl>,
## #   cl_3 <dbl>, cl_3_height <dbl>, precip_12 <dbl>, precip_6 <dbl>,
## #   cl <dbl>, precip <dbl>

A função desc() permite a ordenação decrescente

dados_ordenados <- dados_heathrow %>% 
    arrange(desc(date))

head(dados_ordenados)

## # A tibble: 6 x 26
##   date                usaf  wban  code  station   lat    lon  elev    wd
##   <dttm>              <chr> <chr> <chr> <chr>   <dbl>  <dbl> <dbl> <dbl>
## 1 2019-12-31 23:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25 126. 
## 2 2019-12-31 22:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25 118. 
## 3 2019-12-31 21:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25 110. 
## 4 2019-12-31 20:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25 100. 
## 5 2019-12-31 19:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25 100  
## 6 2019-12-31 18:00:00 0377~ 99999 0377~ HEATHR~  51.5 -0.457    25  93.5
## # ... with 17 more variables: ws <dbl>, ceil_hgt <dbl>, visibility <dbl>,
## #   air_temp <dbl>, dew_point <dbl>, atmos_pres <dbl>, RH <dbl>,
## #   cl_1 <dbl>, cl_1_height <dbl>, cl_2 <dbl>, cl_2_height <dbl>,
## #   cl_3 <dbl>, cl_3_height <dbl>, precip_12 <dbl>, precip_6 <dbl>,
## #   cl <dbl>, precip <dbl>

Funções básicas de tratamento (dplyr) `filter()`:

filter() seleciona colunas de acordo com os seus valores

Funções básicas de tratamento (dplyr) `filter()` (cont.):

Filtrando países (note o operador %in% )

gapminder %>% 
  filter(country %in% c("Brazil", "Argentina", "Chile")) %>% 
  ggplot() +
    geom_line(aes(x = year, y = gdpPercap, color = country )) +
    geom_point(aes(x = year, y = gdpPercap, color = country )) +
    theme_light()

Funções básicas de tratamento (dplyr) `top_n()`:

top_n() seleciona as n linhas maiores de acordo com uma das colunas.

Funções básicas de tratamento (dplyr) `top_n()` (cont.):

Selecionando os países mais ricos em 2007.

Depois aprenderemos como ordenar essas barras

gapminder %>% 
  filter(year == 2007) %>% 
  top_n(5, gdpPercap) %>% 
  ggplot() +
    geom_col(aes(x = country, y = gdpPercap)) +
    theme_light()

Funções básicas de tratamento (dplyr) `top_frac()`:

ricos <- gapminder %>% 
  filter(year == 2007) %>% 
  top_frac(.2, gdpPercap ) %>% 
  mutate(categoria = "Rico")

pobres <-  gapminder %>% 
  filter(year == 2007) %>% 
  top_frac(.2, desc(gdpPercap) ) %>% 
  mutate(categoria = "Pobre")

bind_rows(ricos, pobres) %>% 
  ggplot(aes(x = gdpPercap, y = lifeExp)) +
  geom_point( aes(color = continent)) +
  facet_grid(. ~ categoria, scales = "free_x") +
  geom_smooth(method = "lm", se = FALSE) +
  theme_light()

Funções básicas de tratamento (dplyr) `slice()`:

slice() seleciona “linhas” de um data frame/tibble

Funções básicas de tratamento (dplyr) `slice()`:

classificacao_brasileirao <- read_html("https://pt.wikipedia.org/wiki/Campeonato_Brasileiro_de_Futebol_de_2019_-_S%C3%A9rie_A") %>% 
  html_nodes("table") %>% 
  extract2(7) %>% 
  html_table()


limbo <- classificacao_brasileirao %>% 
  slice(12:16) %>% 
  select(time = Equipes )

limbo

##               time
## 1    Vasco da Gama
## 2 Atlético Mineiro
## 3       Fluminense
## 4         Botafogo
## 5            Ceará

Funções básicas de tratamento (dplyr) `group_by()`:

A função group_by() será bastante usada.

Quem conhece SQL pode estranhar um pouco o comportamento desta função, pois ela não agrupa os dados diminuindo o número de linhas imediatamente.

Mas veja que ela indica que há agrupamento

Ela particiona o tibble. As operações passam a ser executadas em cada partição.

gini_agrupado <- gini %>% 
    select(country, date, value) %>% 
    group_by(country) 
    
gini_agrupado

## # A tibble: 682 x 3
## # Groups:   country [152]
##    country   date  value
##  * <chr>     <chr> <dbl>
##  1 Albania   2012   29  
##  2 Albania   2008   30  
##  3 Algeria   2011   27.6
##  4 Angola    2008   42.7
##  5 Argentina 2017   41.2
##  6 Argentina 2016   42  
##  7 Argentina 2014   41.7
##  8 Argentina 2013   41  
##  9 Argentina 2012   41.4
## 10 Argentina 2011   42.7
## # ... with 672 more rows

Funções básicas de tratamento (dplyr) `group_by()` (cont.):

Para várias operações, o agrupamento faz com que o comportamento seja diferente.

Uma operação bastante usada é numerar as linhas de um tibble.

No tibble agrupado, essa operação acontece em cada grupo.

gini_agrupado %<>% arrange(country, date) %>% 
    mutate(linha = row_number())


datatable(gini_agrupado)

Funções básicas de tratamento (dplyr) `group_by()` (cont.):

As funções lag() no contexto da group_by e lead() no contexto da group_by funcionam dentro de cada grupo (o primeiro value de um grupo não acessa o valor do outro grupo com lag() .

gini_agrupado %<>% mutate(
    value_ant = lag(value),
    delta_value = value - value_ant
    )

datatable(gini_agrupado)

Funções básicas de tratamento (dplyr) `summarise()` (cont.):

Funções básicas de tratamento (dplyr) `group_by()` (cont.):

A função group_by só leva a uma sumarização, ou seja, só transforma o tibble em um tibble com o número de linhas igual ao número de grupos, quando executamos a função summarise()

Se executarmos summarise() sem particionar o tibble, a operação resulta em uma linha.

maiores_temp_dia <- dados_heathrow %>% 
    group_by(date(date)) %>% 
    summarise(
        maxima = max(air_temp),
        minima = min(air_temp),
        media = mean(air_temp)
    )


datatable(maiores_temp_dia) %>% 
    formatRound(c("maxima", "minima", "media"), 1)

A função top_n() no contexto da group_by() retorna os n maiores valores. Se o tibble estiver agrupado, pra cada grupo.

maiores_temp_dia <- dados_heathrow %>% 
    group_by(date(date)) %>% 
    top_n(1, air_temp) %>% 
    ungroup() %>% 
    mutate(
        hora = hour(date),
        estacao = 
            case_when(
                month(date) %in% 1:3 ~ "Inverno",
                month(date) %in% 7:9 ~ "Verão",
                TRUE ~ "Outono/Primavera"
                
            )
    ) %>% 
    select(hora, estacao, air_temp)

ggplot(maiores_temp_dia) +
    geom_density( aes(x = hora, color = estacao )) +
    theme_light()

Leitura de dados de arquivos

Até agora acessamos dados que estavam disponíveis em bibliotecas, mas muitas vezes encontramos dados em arquivos.

De modo geral, as funções da biblioteca readr são mais rápidas do que as da biblioteca base, e também mostram barra de progresso no console. É possível reconhecê-las pelo _ ao invés de .

Leitura de dados de arquivos

O caso mais comum é ler dados em formato de tabela para um tibble

Fonte: (RStudio 2019a)

Leitura de dados de arquivos - Exemplo: CVM

O portal da CVM é uma das minas de ouro de dados

O código abaixo baixa os dados que ainda não estão na nossa base.

Note que vamos usar um exemplo de programação funcional. Semm uso de loop vamos executar uma função para cada linha do data frame/tibble.

existem <- tibble(arquivo = list.files("dados/fundos"))


salva <- function(endereco, arquivo){
    

    print(endereco)
    conteudo <- read_csv2(endereco)
    
    write_csv(conteudo, paste0("dados/fundos/",arquivo))
    
}



baixa_faltantes <- tibble(data_dado = seq(ymd("2017-01-01"), by = "month", ymd(today()-4) )) %>% 
    mutate(data_formato = stamp_date("999912")(data_dado)) %>% 
    mutate(
        endereco = paste0(
            "http://dados.cvm.gov.br/dados/FI/DOC/INF_DIARIO/DADOS/",
            "inf_diario_fi_",
            data_formato,
            ".csv")) %>% 
    mutate(arquivo = paste("inf_diario_fi_",data_formato,".csv")) %>% 
    anti_join(existem, by = c("arquivo" = "arquivo")) %>% 
    select(-data_formato) %>% 
    mutate(data = pmap(.l = list(endereco, arquivo), salva))


le_arquivo <- function(arquivo){
  
    print(arquivo)
      
    conteudo <- read_csv(arquivo, progress = FALSE)

    conteudo
    
}


todos_os_fundos <- tibble(arquivo = list.files("dados/fundos")) %>%
    mutate(arquivo = paste0("dados/fundos/",arquivo )) %>%
    mutate(data = map(arquivo, le_arquivo)) %>% 
    unnest()

## [1] "dados/fundos/inf_diario_fi_ 201701 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201702 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201703 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201704 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201705 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201706 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201707 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201708 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201709 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201710 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201711 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201712 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201801 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201802 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201803 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201804 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201805 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201806 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201807 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201808 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201809 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201810 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201811 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201812 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201901 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201902 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201903 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201904 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201905 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201906 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201907 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201908 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201909 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201910 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201911 .csv"
## [1] "dados/fundos/inf_diario_fi_ 201912 .csv"
## [1] "dados/fundos/inf_diario_fi_ 202001 .csv"
## [1] "dados/fundos/inf_diario_fi_ 202002 .csv"

head(todos_os_fundos)

## # A tibble: 6 x 9
##   arquivo CNPJ_FUNDO DT_COMPTC  VL_TOTAL VL_QUOTA VL_PATRIM_LIQ CAPTC_DIA
##   <chr>   <chr>      <date>        <dbl>    <dbl>         <dbl>     <dbl>
## 1 dados/~ 00.017.02~ 2017-01-02   1.08e8  2.43e13     108099858         0
## 2 dados/~ 00.017.02~ 2017-01-03   1.08e8  2.43e13     108144909         0
## 3 dados/~ 00.017.02~ 2017-01-04   1.08e8  2.43e13     108188649         0
## 4 dados/~ 00.017.02~ 2017-01-05   1.08e8  2.43e13     108234509         0
## 5 dados/~ 00.017.02~ 2017-01-06   1.08e8  2.43e13     108278954         0
## 6 dados/~ 00.017.02~ 2017-01-09   1.08e8  2.43e13     108321610         0
## # ... with 2 more variables: RESG_DIA <dbl>, NR_COTST <dbl>

cadastro_fundos <- read_csv2("http://dados.cvm.gov.br/dados/FIE/CAD/DADOS/inf_cadastral_fie.csv", locale =  locale(encoding = "latin1") )

cotas_verde <- todos_os_fundos %>% 
    filter(CNPJ_FUNDO == "07.455.507/0001-89" )


cotas_verde

## # A tibble: 779 x 9
##    arquivo CNPJ_FUNDO DT_COMPTC  VL_TOTAL VL_QUOTA VL_PATRIM_LIQ CAPTC_DIA
##    <chr>   <chr>      <date>        <dbl>    <dbl>         <dbl>     <dbl>
##  1 dados/~ 07.455.50~ 2017-01-02  1.31e12  8.61e12 1312862006441    8.22e4
##  2 dados/~ 07.455.50~ 2017-01-03  1.31e12  8.62e12 1311291630745    0.    
##  3 dados/~ 07.455.50~ 2017-01-04  1.30e12  8.58e12 1305299654809    2.42e5
##  4 dados/~ 07.455.50~ 2017-01-05  1.29e12  8.55e12 1300213664841    0.    
##  5 dados/~ 07.455.50~ 2017-01-06  1.30e12  8.56e12 1301706718495    3.41e3
##  6 dados/~ 07.455.50~ 2017-01-09  1.30e12  8.55e12 1300644945751    3.27e3
##  7 dados/~ 07.455.50~ 2017-01-10  1.30e12  8.55e12 1300428226690    5.04e3
##  8 dados/~ 07.455.50~ 2017-01-11  1.30e12  8.53e12 1296492684946    7.08e2
##  9 dados/~ 07.455.50~ 2017-01-12  1.32e12  8.60e12 1332935151782    2.50e8
## 10 dados/~ 07.455.50~ 2017-01-13  1.33e12  8.61e12 1333437691636    2.65e3
## # ... with 769 more rows, and 2 more variables: RESG_DIA <dbl>,
## #   NR_COTST <dbl>

Leitura de conteúdo de páginas web

Para efetuar a leitura de páginas WEB, é necessário conhecer como é a estrutura de uma página web.

Uma página web pode ser representada por uma árvore de objetos, também chamada de DOM (Document Object Model).

Esta árvore de objetos é definida pelo conteúdo de linguagem html que existe na página (e pode ser modificado por scripts em javascript e definições de estilo do CSS).

Podem existir objetos de vários tipos em uma página: links, inputs de dados, tabelas, células de tabelas, parágrafos, cabeçalhos etc.

Leitura de conteúdo de páginas WEB (cont.)

Para retirar da página web o conteúdo de que precisamos, temos que analisar como é esta árvore de objetos e que nós desta árvore nos interessam.

Imagine que queremos buscar dados na página de histórico de preços e taxas dos títulos brasileiros

A tecla F12 do Chrome nos permite ver a árvore DOM da página em que estamos navegando.

Leitura de conteúdo de páginas WEB (cont.)

É possível clicar com o botão direito e inspecionar um elemento específico de forma a saber onde ele está na árvore e que tipo de elemento ele é (mesmo que você saiba pouco de html).

O mais importante é saber que uma tag html que define um elemento tem a sintaxe:

<tipo_elemento nome_atributo=valor_attributo>texto do elemento</tipo_elemento>

Mesmo sem saber hmtl, fica claro que queremos esse tal de atributo href dos tais elementos a seja lá o que diabos isso seja (a é um link e href é o destino do link).

Leitura de conteúdo de páginas WEB (cont.)

A biblioteca rvest possibilita a extração destes elementos.

É possível caminhar pela árvore DOM até os nós desejados e atributos que queremos usando html_nodes() e html_attr().

Munidos de uma função que faz download e salva um arquivo, podemos caminhar pelas planilhas e salvá-las

existem <- tibble(arquivo = list.files("dados/titulos"))


salva_planilha <- function(name, endereco){
    
    arquivo <- endereco
    destino <- name

    print(arquivo)
        
    download.file(
        arquivo, 
        paste0("dados/titulos/",destino,".xls" ),
        mode = "wb" ##PRA ARQUIVOS BINÁRIOS,
        )
    
}


links <- read_html("https://sisweb.tesouro.gov.br/apex/f?p=2031:2:0::::") %>% 
    html_nodes("body") %>% 
    html_nodes("a") %>% 
    html_attr("href") %>% 
    enframe(value = "endereco") %>%
    filter(str_detect(endereco, "cosis/sistd/obtem_arquivo/")) %>%
    mutate(destino = paste0(name,".xls")) %>% 
    anti_join(existem, by = c("destino" = "arquivo")) %>%
    select(-destino) %>% 
    mutate(endereco = paste0("https://sisweb.tesouro.gov.br/apex/",endereco) ) %>% 
    mutate(data = map2(name, endereco, salva_planilha ))

Leitura de conteúdo de planilhas Excel

vamos aproveitar as planilhas que baixamos para ver como podemos efetuar a leitura de planilhas Excel.

Agora vamos organizar as planilhas que lemos do site do tesouro.

Vamos usar a biblioteca read_xl.

Precisamos ler todas as sheets de todos os arquivos em um diretório (onde baixamos os arquivos excel do site do tesouro).

A função abaixo lê as sheets de um arquivo.

Aqui fazemos o primeiro uso de uma funçãoda biblioteca stringr, que trata strings, ou seja, cadeias de caracteres.

A função str_glue() possibilita a criação de novas strings a partir de variáveis já existentes de uma forma expressiva.

le_sheets <- function(arquivo){
    
  excel_sheets(str_glue("dados/titulos/{arquivo}")) 
    
}

A função abaixo lê o conteúdo de cada sheet

le_conteudo_sheet <- function(arquivo, sheet){
    
    read_excel(
        str_glue("dados/titulos/{arquivo}"),
        sheet = sheet,
        skip = 2,
        col_names = FALSE,
        col_types = "text"
    ) 
    
}

Leitura de conteúdo de planilhas Excel (cont.)

Munidos destas duas funções, podemos guardar tudo em um só data frame.

Primeiro, definindo as sheets a ler

sheets_pra_ler <- list.files("dados/titulos") %>% 
    enframe(value = "arquivo") %>%
    mutate(sheet = map(arquivo, le_sheets)) %>% 
    unnest(cols = sheet)

Depois lendo as sheets para um data frame

Neste momento, vamos executar a função em processamento paralelo, usando todos os processadores da nossa máquina.

Para isso, vamos lançar mão da biblioteca furrr%.

Ela contém adaptações de todas as funções map. São as funções baseadas na future_map.

Antes de rodar a future_map, determinamos a forma de paralelização. No nosso caso, plan(multiprocess).

É legal notar a forma como usamos expressão regular aqui

plan(multiprocess)

sheets_lidas <- sheets_pra_ler %>% 
  mutate(titulo = str_extract(sheet,"[^[0-9]]*" )) %>% 
  mutate(vencimento = str_extract(sheet,"[0-9]{6}" )) %>% 
  mutate(vencimento = dmy(vencimento)) %>%
  mutate(
      arquivo_out = arquivo,
      sheet_out = sheet
  ) %>% 
  mutate(data = future_map2(.x = arquivo, .y = sheet , le_conteudo_sheet, .progress = TRUE)) %>% 
  unnest(data)

## 
 Progress: -------------------------------------------------------------------------------  100%
 Progress: -------------------------------------------------------------------------------- 100%

Leitura de conteúdo de planilhas Excel (cont.)

Uma última arrumada

taxas_titulos <- sheets_lidas %>% 
    rename(
        data = 8,
        taxa_compra_manha = 9,
        taxa_venda_manha = 10,
        pu_compra_manha = 11,
        pu_venda_manha = 12,
        pu_base_manha = 13
    ) %>% 
    mutate(
        data = if_else(
            str_detect(data, "/"),
            dmy(data),
            as.numeric(data) + ymd("1899-12-31")
        )
    ) %>% 
    mutate(
        titulo = str_trim(titulo), 
        taxa_compra_manha = as.numeric(taxa_compra_manha),
        taxa_venda_manha = as.numeric(taxa_venda_manha),
        pu_compra_manha = as.numeric(pu_compra_manha),
        pu_venda_manha = as.numeric(pu_venda_manha),
        pu_base_manha = as.numeric(pu_base_manha)
    ) %>% 
    select(
        titulo,
        vencimento,
        data,
        taxa_compra_manha,
        taxa_venda_manha,
        pu_compra_manha,
        pu_venda_manha,
        pu_base_manha
    ) %>% 
    mutate(
        titulo = if_else(
            titulo == "NTN-B Principal", 
            "NTN-B Princ", 
            titulo
        )
    )

Leitura de conteúdo de planilhas Excel (cont.)

Aí fica fácil fazer a análise que desejarmos

ntnb_2045 <- taxas_titulos %>% 
    filter(
        titulo == "NTN-B Princ",
        vencimento == ymd("2045-05-15")
    ) %>% 
    mutate(taxa = (taxa_compra_manha +taxa_venda_manha)/2 )



ggplot(ntnb_2045) +
    geom_line(aes(x = data, y = taxa) ) +
    scale_x_date(
        date_breaks = "3 months",
        limits = c(ymd("2016-12-01"),NA),
        labels = date_format("%m/%y")
    ) +
    scale_y_continuous(
        labels = percent_format(accuracy = 0.1)
    ) + 
    labs(y = "NTN-B Principal 2045", x = "Data") +
    theme_light()

Outro exemplo de leitura de página

Neste caso, precisamos entender quando requisição é feita pela página e emular essas requisições passando os parâmetros que devem ser passados via método POST.

Primeiro usamos a função crossing para criar

url <- "http://www.ceagesp.gov.br/entrepostos/servicos/cotacoes/#cotacao"


le_pagina_ceagesp <- function(grupo, data){

  print(str_glue("Grupo:{grupo}, data:{data}"))

  dados <- NA 
  fd <- list(  
    cot_grupo  = grupo,
    cot_data = data
  )
  
  resp <- POST(url, body=fd, encode="form")
  
  tabela <- content(resp) %>% 
    html_nodes("table") %>% 
    extract2(1) %>% 
    html_table()
  
  nomes <- tabela %>% slice(2)
  
  dados <- tabela %>% slice(3:nrow(tabela)) %>% 
    mutate(data = dmy(data))
  
  names(dados) <- c(nomes, "data_precos")

  dados
  
}


le_pagina_ceagesp_ignora_erro <- possibly(le_pagina_ceagesp, otherwise = NA) 


grupos <- c(
  "FRUTAS",
  "LEGUMES",
  "VERDURAS",
  "DIVERSOS",
  "FLORES",
  "PESCADOS"
)

datas <- seq(from = today()-3, by = 1, to = today()-1) %>% 
  format("%d/%m/%Y") 

dados_ceagesp <-  enframe(grupos, value = "grupo", name = "id_grupo") %>%
  crossing(enframe(datas, value = "data", name = "id_data")) %>% 
  mutate(leitura = map2(.x = grupo, .y = data, le_pagina_ceagesp_ignora_erro )) %>% 
  filter(!is.na(leitura)) %>%
  unnest(cols = leitura)

## Grupo:FRUTAS, data:28/02/2020
## Grupo:FRUTAS, data:29/02/2020
## Grupo:FRUTAS, data:01/03/2020
## Grupo:LEGUMES, data:28/02/2020
## Grupo:LEGUMES, data:29/02/2020
## Grupo:LEGUMES, data:01/03/2020
## Grupo:VERDURAS, data:28/02/2020
## Grupo:VERDURAS, data:29/02/2020
## Grupo:VERDURAS, data:01/03/2020
## Grupo:DIVERSOS, data:28/02/2020
## Grupo:DIVERSOS, data:29/02/2020
## Grupo:DIVERSOS, data:01/03/2020
## Grupo:FLORES, data:28/02/2020
## Grupo:FLORES, data:29/02/2020
## Grupo:FLORES, data:01/03/2020
## Grupo:PESCADOS, data:28/02/2020
## Grupo:PESCADOS, data:29/02/2020
## Grupo:PESCADOS, data:01/03/2020

head(dados_ceagesp)

## # A tibble: 6 x 12
##   id_grupo grupo id_data data  Produto Classificação `Uni/Peso` Menor Comun
##      <int> <chr>   <int> <chr> <chr>   <chr>         <chr>      <chr> <chr>
## 1        1 FRUT~       1 28/0~ ABACAT~ A             KG         12,28 14,08
## 2        1 FRUT~       1 28/0~ ABACAT~ A BOCA 8 a 10 KG         2,36  2,53 
## 3        1 FRUT~       1 28/0~ ABACAT~ B BOCA 11 e ~ KG         1,9   2,08 
## 4        1 FRUT~       1 28/0~ ABACAT~ A BOCA 8 a 10 KG         1,62  1,81 
## 5        1 FRUT~       1 28/0~ ABACAT~ B BOCA 11 e ~ KG         1,23  1,42 
## 6        1 FRUT~       1 28/0~ ABACAT~ A BOCA 8 a 10 KG         1,6   1,79 
## # ... with 3 more variables: Maior <chr>, Quilo <chr>, data_precos <date>

Leitura de página como um robô

Em dois casos precisamos emular um browser para obter o conteúdo html que desejamos, para então usar a rvest

Quando a página precisa executar um script javascript para que fique pronta
Quando a página exige muitos passos de navegação antes de chegarmos aos dados desejados;

Usando Selenium

Selenium é um famoso automatizador de browser. Através dele, é possível enviar cliques, preencher caixas de texto etc. até se chegar aos dados desejados. Além disso, ele executa os scripts javascripts que a página exige.

A biblioteca RSelenium possibilita o uso na linguagem R

Abaixo um exemplo de ativação do firefox a partir do R com RSelenium

rs <- rsDriver(
    browser = "firefox",
    extraCapabilities = list(
        `mox:firefoxOptions` = list(
            binary = "C:/Program Files (x86)/Mozilla Firefox/firefox.exe"
        )
    )
)

Baixando dados da ANEEL

library(RSelenium)
library(rvest)
library(httr)
library(magrittr)

# initiate selenium     -- start it this way every time
# make sure the location of firefox on your machine is correct
rs <- rsDriver(
  browser = "firefox", 
  extraCapabilities = list(
    `mox:firefoxOptions` = list(
      binary = "C:/Program Files/Mozilla Firefox/firefox.exe"
    )
  )
)


preparou <-  function(){
  outputzipfile <- "nada"
  try(outputzipfile <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[2]/div[2]/div/div[1]"))  
  if(typeof(outputzipfile) == "S4"){
    resposta <- outputzipfile$getElementText() == "Output Zip File"
  }
  else{
    resposta <-  FALSE
  }
  
  resposta
  
}


rsc <- rs$client

rsc$navigate("https://sigel.aneel.gov.br/Down/")

imagem <- rsc$findElement(using = "css selector",  "#uniqName_9_1 > div:nth-child(1)" )

imagem$clickElement()

Sys.sleep(3)

UHE <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[1]/div[1]/div[1]/div[2]/div/div[3]/div[1]")

UHE$clickElement()

GD <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[1]/div[1]/div[1]/div[2]/div/div[1]/div[1]")

GD$clickElement()

baixar <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[1]/div[2]/div[2]/div")

baixar$clickElement()

while(!preparou()){
  
  print("esperando")
  Sys.sleep(1)
  
}
  
arquivo_download <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[2]/div[2]/div/div[2]/div/a")

endereco <- arquivo_download$getElementText()

download.file(url = unlist(endereco), destfile = "c:\\temp\\gd.zip")

E quando o desenvolvedor da página web tá manguaçado?

Nem todas as páginas são fáceis de ler.

A página que mostra os DIs na BMF, por exemplo é esquisita:

O nosso objeto de interesse aparece na ferramenta chamada por F12, mas não é encontrada pela rvest ao ler o HTML

Lendo páginas difíceis

Existem duas técnicas para o Web Scraping:

Simular a navegação num browser emulado (fora do escopo deste curso)
Interpretar páginas com base no conteúdo recebido

Na página que visualizamos da BMF, o conteúdo recebido inclui um script javascript que popula a tabela, por isso não conseguimos reconhecer o conteúdo da tabela de pronto, pois ela só é carregada pela execução do script.

Lendo páginas difíceis (cont.)

Conseguimos, no entanto, extrair do script as informações de que precisamos

Lendo páginas difíceis (cont.)

O código abaixo extrai da página as informações de que precisamos.

Para extrair as informações, ele usa expressões regulares.

Repare na função que silencia um possível erro. Este erro pode acontecer se passarmos um dia que não tem dados. Precisamos disso pois vamos buscar todas as datas possíveis.

Essa não é a melhor forma de tratar um erro deste tipo. Veremos mais adiante

cotacoes <- tibble(
    id = 0:10,
    cotacao = c(
        "ajuste_ant",
        "ajuste_corrig",
        "preco_abert",
        "preco_min",
        "preco_max",
        "preco_med",
        "ult_preco",
        "ajuste",
        "var_pontos",
        "ult_of_compra",
        "ult_of_venda"
    )
)

#QUEM INVENTOU ISSO?
siglas_mes <- tibble(
    sigla_mes = c(
        "F",
        "G",
        "H",
        "J",
        "K",
        "M",
        "N",
        "Q",
        "U",
        "V",
        "X",
        "Z"
        ),
    mes = 1:12
)



pagina <- NA

try( 
    {pagina <- read_html("http://www2.bmf.com.br/pages/portal/bmfbovespa/boletim1/SistemaPregao1.asp?pagetype=pop&caminho=Resumo%20Estat%EDstico%20-%20Sistema%20Preg%E3o&Data=10/01/2017&Mercadoria=DI1") %>% 
    html_nodes("script") %>% 
    extract2(13) %>% 
    html_text()}
    ,
    silent = TRUE
)


linhas_impares <- str_extract_all(pagina, "MercFut1 \\+ '<td ALIGN=\"right\" CLASS=\"tabelaConteudo1\">[0-9.,]*") %>% 
    extract2(1) %>% 
    enframe() %>% 
    mutate(
        ativo = (row_number()-1) %/% 11 * 2,
        tipo_cotacao = (row_number() - 1) %% 11
    )


linhas_pares <- str_extract_all(pagina, "MercFut1 \\+ '<td ALIGN=\"right\" CLASS=\"tabelaConteudo2\">[0-9.,]*") %>% 
    extract2(1) %>% 
    enframe() %>% 
    mutate(
        ativo = (row_number()-1) %/% 11 * 2 + 1,
        tipo_cotacao = (row_number() - 1) %% 11
    )


linhas <- linhas_impares %>% 
    bind_rows(linhas_pares) %>%
    arrange(ativo) %>% 
    mutate(valor = str_extract_all(value, ">[0-9.,]*")) %>% 
    mutate(valor = str_sub(valor, 2)) %>% 
    select(ativo, tipo_cotacao, valor)



ativos <- str_extract_all(pagina, "MercFut3 = MercFut3 \\+ '</tr><td ALIGN=\"center\" CLASS=\"tabelaConteudo[0-9]\">[A-Z][0-9][0-9]") %>%
    extract2(1) %>% 
    enframe() %>% 
    mutate(
        nome_ativo = str_sub(value,-3),
        id = name - 1
    ) %>% 
    select(
        id,
        nome_ativo
    )


linhas %<>% left_join(ativos, by = c("ativo" = "id")) %>% 
    select(nome_ativo, tipo_cotacao, valor) %>% 
    mutate(
        sigla_mes = str_sub(nome_ativo, 1,1),
        ano = as.numeric(str_sub(nome_ativo, 2,4)) + 2000
    ) %>% 
    left_join(cotacoes, by = c("tipo_cotacao" = "id") ) %>% 
    left_join(siglas_mes, by = c("sigla_mes")) %>% 
    mutate(vencimento = make_date(ano, mes, 1)) %>%
    mutate(
        valor = parse_number(
            valor, 
            locale = locale(
                decimal_mark = ",", 
                grouping_mark = "." 
            )
        )
    ) %>% 
    select(
        nome_ativo,
        vencimento,
        cotacao,
        valor
    )

datatable(linhas)

Lendo páginas difíceis (cont.)

Agora vamos executar para todos os dias desde janeiro de 2017.

Preparamos a função…

le_uma_pagina_bmf <- function(dados){
    
    
    
    data <- pull(dados, data_in) %>% 
        stamp_date("31/01/2017")(.)
    

    
    pagina <- NA

    try( 
        {pagina <- read_html(paste0("http://www2.bmf.com.br/pages/portal/bmfbovespa/boletim1/SistemaPregao1.asp?pagetype=pop&caminho=Resumo%20Estat%EDstico%20-%20Sistema%20Preg%E3o&Data=",data,"&Mercadoria=DI1")) %>% 
        html_nodes("script") %>% 
        extract2(13) %>% 
        html_text()}
        ,
        silent = TRUE
    )
        
    
        
    if (!is.na(pagina)){
        
        
        
        linhas_impares <- str_extract_all(pagina, "MercFut1 \\+ '<td ALIGN=\"right\" CLASS=\"tabelaConteudo1\">[0-9.,]*") %>% 
            extract2(1) %>% 
            enframe() %>% 
            mutate(
                ativo = (row_number()-1) %/% 11 * 2,
                tipo_cotacao = (row_number() - 1) %% 11
            )
        
        
        linhas_pares <- str_extract_all(pagina, "MercFut1 \\+ '<td ALIGN=\"right\" CLASS=\"tabelaConteudo2\">[0-9.,]*") %>% 
            extract2(1) %>% 
            enframe() %>% 
            mutate(
                ativo = (row_number()-1) %/% 11 * 2 + 1,
                tipo_cotacao = (row_number() - 1) %% 11
            )
        
        
        linhas <- linhas_impares %>% 
            bind_rows(linhas_pares) %>%
            arrange(ativo) %>% 
            mutate(valor = str_extract_all(value, ">[0-9.,]*")) %>% 
            mutate(valor = str_sub(valor, 2)) %>% 
            select(ativo, tipo_cotacao, valor)
        
        
        
        ativos <- str_extract_all(pagina, "MercFut3 = MercFut3 \\+ '</tr><td ALIGN=\"center\" CLASS=\"tabelaConteudo[0-9]\">[A-Z][0-9][0-9]") %>%
            extract2(1) %>% 
            enframe() %>% 
            mutate(
                nome_ativo = str_sub(value,-3),
                id = name - 1
            ) %>% 
            select(
                id,
                nome_ativo
            )
        
        
        linhas %<>% left_join(ativos, by = c("ativo" = "id")) %>% 
            select(nome_ativo, tipo_cotacao, valor) %>% 
            mutate(
                sigla_mes = str_sub(nome_ativo, 1,1),
                ano = as.numeric(str_sub(nome_ativo, 2,4)) + 2000
            ) %>% 
            left_join(cotacoes, by = c("tipo_cotacao" = "id") ) %>% 
            left_join(siglas_mes, by = c("sigla_mes")) %>% 
            mutate(vencimento = make_date(ano, mes, 1)) %>%
            mutate(
                valor = parse_number(
                    valor, 
                    locale = locale(
                        decimal_mark = ",", 
                        grouping_mark = "." 
                    )
                )
            ) %>% 
            select(
                nome_ativo,
                vencimento,
                cotacao,
                valor
            )
        
        linhas
    }
    else
    {
        tibble(dia_inutil = TRUE )
    }
}

E executamos para todos os dias.

dados_todas_as_datas <- seq.Date(
        ymd("2020-01-01"), 
        to = today(), 
        by = "day") %>% 
    enframe() %>% 
    rename(data_out = value) %>% 
    mutate(data_in = data_out) %>% 
    group_by(data_out, name) %>% 
    nest_legacy() %>% 
    mutate(data = map(data, le_uma_pagina_bmf)) %>% 
    unnest_legacy()

Lendo páginas difíceis (cont.)

O data frame com todos os dados ficou assim.

Veja como as funções da biblioteca kable e kableExtra dão mais controle na criação das tabelas.

Note como ele ainda não está “Tidy”. Por quê?

dados_todas_as_datas %>%
    ungroup() %>% 
    filter(is.na(dia_inutil)) %>% 
    select(-dia_inutil, - name) %>% 
    head(n = 200) %>% 
    mutate(
        data_out = stamp_date("31/12/2010")(data_out),
        vencimento = stamp_date("12/%Y")(vencimento),
        valor = number(valor, accuracy = 0.001, decimal.mark = ",", big.mark = ".")
    ) %>% 
    kable(
        col.names = 
            c(
                "Data",
                "Ativo",
                "Vencimento",
                "Tipo Cotação",
                "Cotação"
            ),
        align = 
            c("l","l","l","l","r")
    ) %>% 
    kable_styling(
        
    ) %>% 
    row_spec(
        seq(2, 200, 2),
        background = "#eeeeee"
    )

Data	Ativo	Vencimento	Tipo Cotação	Cotação
02/01/2020	F20	01/2020	ajuste_ant	99.999,960
02/01/2020	F20	01/2020	ajuste_corrig	99.999,960
02/01/2020	F20	01/2020	preco_abert	0,000
02/01/2020	F20	01/2020	preco_min	0,000
02/01/2020	F20	01/2020	preco_max	0,000
02/01/2020	F20	01/2020	preco_med	0,000
02/01/2020	F20	01/2020	ult_preco	0,000
02/01/2020	F20	01/2020	ajuste	100.000,000
02/01/2020	F20	01/2020	var_pontos	0,040
02/01/2020	F20	01/2020	ult_of_compra	0,000
02/01/2020	F20	01/2020	ult_of_venda	0,000
02/01/2020	G20	02/2020	ajuste_ant	99.623,950
02/01/2020	G20	02/2020	ajuste_corrig	99.623,950
02/01/2020	G20	02/2020	preco_abert	4,410
02/01/2020	G20	02/2020	preco_min	4,409
02/01/2020	G20	02/2020	preco_max	4,422
02/01/2020	G20	02/2020	preco_med	4,410
02/01/2020	G20	02/2020	ult_preco	4,410
02/01/2020	G20	02/2020	ajuste	99.623,790
02/01/2020	G20	02/2020	var_pontos	0,160
02/01/2020	G20	02/2020	ult_of_compra	4,409
02/01/2020	G20	02/2020	ult_of_venda	4,410
02/01/2020	H20	03/2020	ajuste_ant	99.322,160
02/01/2020	H20	03/2020	ajuste_corrig	99.322,160
02/01/2020	H20	03/2020	preco_abert	4,370
02/01/2020	H20	03/2020	preco_min	4,370
02/01/2020	H20	03/2020	preco_max	4,380
02/01/2020	H20	03/2020	preco_med	4,378
02/01/2020	H20	03/2020	ult_preco	4,375
02/01/2020	H20	03/2020	ajuste	99.322,620
02/01/2020	H20	03/2020	var_pontos	0,460
02/01/2020	H20	03/2020	ult_of_compra	4,365
02/01/2020	H20	03/2020	ult_of_venda	4,375
02/01/2020	J20	04/2020	ajuste_ant	98.959,420
02/01/2020	J20	04/2020	ajuste_corrig	98.959,420
02/01/2020	J20	04/2020	preco_abert	4,320
02/01/2020	J20	04/2020	preco_min	4,320
02/01/2020	J20	04/2020	preco_max	4,338
02/01/2020	J20	04/2020	preco_med	4,335
02/01/2020	J20	04/2020	ult_preco	4,336
02/01/2020	J20	04/2020	ajuste	98.961,120
02/01/2020	J20	04/2020	var_pontos	1,700
02/01/2020	J20	04/2020	ult_of_compra	4,338
02/01/2020	J20	04/2020	ult_of_venda	4,340
02/01/2020	K20	05/2020	ajuste_ant	98.628,140
02/01/2020	K20	05/2020	ajuste_corrig	98.628,140
02/01/2020	K20	05/2020	preco_abert	4,350
02/01/2020	K20	05/2020	preco_min	4,305
02/01/2020	K20	05/2020	preco_max	4,350
02/01/2020	K20	05/2020	preco_med	4,312
02/01/2020	K20	05/2020	ult_preco	4,310
02/01/2020	K20	05/2020	ajuste	98.636,310
02/01/2020	K20	05/2020	var_pontos	8,170
02/01/2020	K20	05/2020	ult_of_compra	4,300
02/01/2020	K20	05/2020	ult_of_venda	4,350
02/01/2020	M20	06/2020	ajuste_ant	98.296,680
02/01/2020	M20	06/2020	ajuste_corrig	98.296,680
02/01/2020	M20	06/2020	preco_abert	4,300
02/01/2020	M20	06/2020	preco_min	4,300
02/01/2020	M20	06/2020	preco_max	4,305
02/01/2020	M20	06/2020	preco_med	4,300
02/01/2020	M20	06/2020	ult_preco	4,305
02/01/2020	M20	06/2020	ajuste	98.304,620
02/01/2020	M20	06/2020	var_pontos	7,940
02/01/2020	M20	06/2020	ult_of_compra	4,305
02/01/2020	M20	06/2020	ult_of_venda	4,340
02/01/2020	N20	07/2020	ajuste_ant	97.952,780
02/01/2020	N20	07/2020	ajuste_corrig	97.952,780
02/01/2020	N20	07/2020	preco_abert	4,320
02/01/2020	N20	07/2020	preco_min	4,315
02/01/2020	N20	07/2020	preco_max	4,325
02/01/2020	N20	07/2020	preco_med	4,318
02/01/2020	N20	07/2020	ult_preco	4,315
02/01/2020	N20	07/2020	ajuste	97.959,160
02/01/2020	N20	07/2020	var_pontos	6,380
02/01/2020	N20	07/2020	ult_of_compra	4,315
02/01/2020	N20	07/2020	ult_of_venda	4,320
02/01/2020	Q20	08/2020	ajuste_ant	97.556,440
02/01/2020	Q20	08/2020	ajuste_corrig	97.556,440
02/01/2020	Q20	08/2020	preco_abert	4,350
02/01/2020	Q20	08/2020	preco_min	4,345
02/01/2020	Q20	08/2020	preco_max	4,350
02/01/2020	Q20	08/2020	preco_med	4,346
02/01/2020	Q20	08/2020	ult_preco	4,345
02/01/2020	Q20	08/2020	ajuste	97.568,640
02/01/2020	Q20	08/2020	var_pontos	12,200
02/01/2020	Q20	08/2020	ult_of_compra	4,315
02/01/2020	Q20	08/2020	ult_of_venda	4,355
02/01/2020	U20	09/2020	ajuste_ant	97.205,500
02/01/2020	U20	09/2020	ajuste_corrig	97.205,500
02/01/2020	U20	09/2020	preco_abert	0,000
02/01/2020	U20	09/2020	preco_min	0,000
02/01/2020	U20	09/2020	preco_max	0,000
02/01/2020	U20	09/2020	preco_med	0,000
02/01/2020	U20	09/2020	ult_preco	0,000
02/01/2020	U20	09/2020	ajuste	97.221,350
02/01/2020	U20	09/2020	var_pontos	15,850
02/01/2020	U20	09/2020	ult_of_compra	4,340
02/01/2020	U20	09/2020	ult_of_venda	4,370
02/01/2020	V20	10/2020	ajuste_ant	96.828,170
02/01/2020	V20	10/2020	ajuste_corrig	96.828,170
02/01/2020	V20	10/2020	preco_abert	4,405
02/01/2020	V20	10/2020	preco_min	4,380
02/01/2020	V20	10/2020	preco_max	4,410
02/01/2020	V20	10/2020	preco_med	4,393
02/01/2020	V20	10/2020	ult_preco	4,380
02/01/2020	V20	10/2020	ajuste	96.849,060
02/01/2020	V20	10/2020	var_pontos	20,890
02/01/2020	V20	10/2020	ult_of_compra	4,380
02/01/2020	V20	10/2020	ult_of_venda	4,390
02/01/2020	X20	11/2020	ajuste_ant	96.454,660
02/01/2020	X20	11/2020	ajuste_corrig	96.454,660
02/01/2020	X20	11/2020	preco_abert	0,000
02/01/2020	X20	11/2020	preco_min	0,000
02/01/2020	X20	11/2020	preco_max	0,000
02/01/2020	X20	11/2020	preco_med	0,000
02/01/2020	X20	11/2020	ult_preco	0,000
02/01/2020	X20	11/2020	ajuste	96.479,550
02/01/2020	X20	11/2020	var_pontos	24,890
02/01/2020	X20	11/2020	ult_of_compra	4,415
02/01/2020	X20	11/2020	ult_of_venda	4,445
02/01/2020	Z20	12/2020	ajuste_ant	96.074,030
02/01/2020	Z20	12/2020	ajuste_corrig	96.074,030
02/01/2020	Z20	12/2020	preco_abert	0,000
02/01/2020	Z20	12/2020	preco_min	0,000
02/01/2020	Z20	12/2020	preco_max	0,000
02/01/2020	Z20	12/2020	preco_med	0,000
02/01/2020	Z20	12/2020	ult_preco	0,000
02/01/2020	Z20	12/2020	ajuste	96.102,390
02/01/2020	Z20	12/2020	var_pontos	28,360
02/01/2020	Z20	12/2020	ult_of_compra	4,465
02/01/2020	Z20	12/2020	ult_of_venda	4,490
02/01/2020	F21	01/2021	ajuste_ant	95.654,620
02/01/2020	F21	01/2021	ajuste_corrig	95.654,620
02/01/2020	F21	01/2021	preco_abert	4,530
02/01/2020	F21	01/2021	preco_min	4,520
02/01/2020	F21	01/2021	preco_max	4,560
02/01/2020	F21	01/2021	preco_med	4,533
02/01/2020	F21	01/2021	ult_preco	4,520
02/01/2020	F21	01/2021	ajuste	95.687,700
02/01/2020	F21	01/2021	var_pontos	33,080
02/01/2020	F21	01/2021	ult_of_compra	4,520
02/01/2020	F21	01/2021	ult_of_venda	4,525
02/01/2020	J21	04/2021	ajuste_ant	94.447,570
02/01/2020	J21	04/2021	ajuste_corrig	94.447,570
02/01/2020	J21	04/2021	preco_abert	4,700
02/01/2020	J21	04/2021	preco_min	4,680
02/01/2020	J21	04/2021	preco_max	4,730
02/01/2020	J21	04/2021	preco_med	4,695
02/01/2020	J21	04/2021	ult_preco	4,690
02/01/2020	J21	04/2021	ajuste	94.483,390
02/01/2020	J21	04/2021	var_pontos	35,820
02/01/2020	J21	04/2021	ult_of_compra	4,690
02/01/2020	J21	04/2021	ult_of_venda	4,700
02/01/2020	N21	07/2021	ajuste_ant	93.129,700
02/01/2020	N21	07/2021	ajuste_corrig	93.129,700
02/01/2020	N21	07/2021	preco_abert	4,890
02/01/2020	N21	07/2021	preco_min	4,880
02/01/2020	N21	07/2021	preco_max	4,930
02/01/2020	N21	07/2021	preco_med	4,899
02/01/2020	N21	07/2021	ult_preco	4,880
02/01/2020	N21	07/2021	ajuste	93.172,860
02/01/2020	N21	07/2021	var_pontos	43,160
02/01/2020	N21	07/2021	ult_of_compra	4,870
02/01/2020	N21	07/2021	ult_of_venda	4,890
02/01/2020	V21	10/2021	ajuste_ant	91.679,300
02/01/2020	V21	10/2021	ajuste_corrig	91.679,300
02/01/2020	V21	10/2021	preco_abert	5,120
02/01/2020	V21	10/2021	preco_min	5,080
02/01/2020	V21	10/2021	preco_max	5,120
02/01/2020	V21	10/2021	preco_med	5,091
02/01/2020	V21	10/2021	ult_preco	5,080
02/01/2020	V21	10/2021	ajuste	91.729,850
02/01/2020	V21	10/2021	var_pontos	50,550
02/01/2020	V21	10/2021	ult_of_compra	5,070
02/01/2020	V21	10/2021	ult_of_venda	5,090
02/01/2020	F22	01/2022	ajuste_ant	90.251,950
02/01/2020	F22	01/2022	ajuste_corrig	90.251,950
02/01/2020	F22	01/2022	preco_abert	5,240
02/01/2020	F22	01/2022	preco_min	5,230
02/01/2020	F22	01/2022	preco_max	5,300
02/01/2020	F22	01/2022	preco_med	5,274
02/01/2020	F22	01/2022	ult_preco	5,250
02/01/2020	F22	01/2022	ajuste	90.292,140
02/01/2020	F22	01/2022	var_pontos	40,190
02/01/2020	F22	01/2022	ult_of_compra	5,240
02/01/2020	F22	01/2022	ult_of_venda	5,250
02/01/2020	J22	04/2022	ajuste_ant	88.832,120
02/01/2020	J22	04/2022	ajuste_corrig	88.832,120
02/01/2020	J22	04/2022	preco_abert	5,440
02/01/2020	J22	04/2022	preco_min	5,410
02/01/2020	J22	04/2022	preco_max	5,460
02/01/2020	J22	04/2022	preco_med	5,440
02/01/2020	J22	04/2022	ult_preco	5,410
02/01/2020	J22	04/2022	ajuste	88.876,780
02/01/2020	J22	04/2022	var_pontos	44,660
02/01/2020	J22	04/2022	ult_of_compra	5,400
02/01/2020	J22	04/2022	ult_of_venda	5,410
02/01/2020	N22	07/2022	ajuste_ant	87.456,380
02/01/2020	N22	07/2022	ajuste_corrig	87.456,380

Pivoteando

Reparamos que o data frame no slide anterior não está “Tidy”, ou seja não está de forma que cada linha represente um evento e cada coluna represente um atributo do evento.

Para nós, neste caso, um evento é formado por todas as informações de um ativo em um dia e não uma só das informações de um ativo em um dia.

Isso porque é extremamente comum fazermos contas com mais de uma informação do ativo em um dia (máxima - mínima, por exemplo)

Tidyr

O pacote Tidyr ajuda a arrumar os data frames dessas formas. O hex sticker dele é bem explicativo.

Tidyr - `pivot_longer()` e `pivot_wider()`

Os nomes das principais funções mudaram em setembro de 2019 (quando saiu a versão 1.0.0). Antes se chamavam gather() e spread() e agora se chamam pivot_wider() e pivot_longer(), o que é mais intuitivo.

Pivoteando nosso data frame com muitas linhas

O nosso data frame tem cada tipo de cotação em cada linha e gostaríamos que essas linhas fossem transformadas em colunas.

A função que faz isso se chama pivot_wider()

Seus parâmetros mais usados são:

data, que é o data frame a ser tratado
names_from, que é o atributo de onde vêm os nomes para os novos atributos
values_from, o atributo de onde vêm os valores dos novos atributos

dados_todas_as_datas %>% 
    pivot_wider(
        names_from = cotacao,
        values_from = valor
    ) %>% 
    str()

## Classes 'tbl_df', 'tbl' and 'data.frame':    1555 obs. of  17 variables:
##  $ data_out     : Date, format: "2020-01-01" "2020-01-02" ...
##  $ name         : int  1 2 2 2 2 2 2 2 2 2 ...
##  $ dia_inutil   : logi  TRUE NA NA NA NA NA ...
##  $ nome_ativo   : chr  NA "F20" "G20" "H20" ...
##  $ vencimento   : Date, format: NA "2020-01-01" ...
##  $ NA           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ ajuste_ant   : num  NA 100000 99624 99322 98959 ...
##  $ ajuste_corrig: num  NA 100000 99624 99322 98959 ...
##  $ preco_abert  : num  NA 0 4.41 4.37 4.32 4.35 4.3 4.32 4.35 0 ...
##  $ preco_min    : num  NA 0 4.41 4.37 4.32 ...
##  $ preco_max    : num  NA 0 4.42 4.38 4.34 ...
##  $ preco_med    : num  NA 0 4.41 4.38 4.33 ...
##  $ ult_preco    : num  NA 0 4.41 4.38 4.34 ...
##  $ ajuste       : num  NA 100000 99624 99323 98961 ...
##  $ var_pontos   : num  NA 0.04 0.16 0.46 1.7 ...
##  $ ult_of_compra: num  NA 0 4.41 4.36 4.34 ...
##  $ ult_of_venda : num  NA 0 4.41 4.38 4.34 ...

Pivoteando nosso data frame com muitas colunas

É muito comum recebermos os dados com colunas que deviam ser valores de um atributo, e não um atributo em si.

O exemplo clássico é a colocação de datas nas colunas do dado, como nos dados retirados do site Datasus

read_csv2(
    "dados/siab/cadastro_numero_familias.csv", 
    skip = 3,
    locale = locale(encoding = "latin1" ) 
    ) %>% 
    glimpse()

## Observations: 5,502
## Variables: 19
## $ Município <chr> "110001 Alta Floresta D'Oeste", "110037 Alto Alegre ...
## $ `1998`    <chr> "2", "1224", "-", "3797", "676", "-", "-", "-", "529...
## $ `1999`    <chr> "3458", "2204", "546", "2763", "9922", "2268", "-", ...
## $ `2000`    <chr> "4668", "1911", "2002", "2832", "11790", "2615", "17...
## $ `2001`    <chr> "5078", "1994", "2348", "3924", "11648", "3061", "18...
## $ `2002`    <chr> "5078", "1976", "2348", "3907", "11654", "3732", "18...
## $ `2003`    <chr> "5077", "1804", "2021", "3463", "9931", "5499", "138...
## $ `2004`    <chr> "4925", "1766", "932", "4077", "18830", "7580", "180...
## $ `2005`    <chr> "5865", "1698", "1956", "4141", "14531", "7470", "18...
## $ `2006`    <chr> "6518", "3623", "1994", "4244", "16509", "6459", "19...
## $ `2007`    <chr> "6489", "3500", "2922", "1642", "17094", "7502", "19...
## $ `2008`    <chr> "7093", "3682", "3195", "4493", "16936", "7712", "23...
## $ `2009`    <chr> "6946", "3172", "3210", "4693", "16948", "7546", "23...
## $ `2010`    <chr> "6836", "3292", "1480", "1586", "16528", "7112", "25...
## $ `2011`    <chr> "6918", "3416", "2551", "3122", "16912", "7378", "18...
## $ `2012`    <chr> "6699", "3448", "3175", "4162", "18109", "7385", "19...
## $ `2013`    <chr> "6336", "3617", "3029", "4242", "19270", "9730", "19...
## $ `2014`    <chr> "6201", "3343", "3438", "4242", "14661", "9084", "19...
## $ `2015`    <chr> "6181", "3280", "-", "-", "6307", "-", "-", "-", "19...

Pivoteando nosso data frame com muitas colunas (cont.)

Acontece que “2009” não é um atributo, mas sim o valor de um atributo que deveria ser data

A função pivot_longer() faz a operação de que precisamos.

data, que é o data frame a ser tratado
names_from, que é o atributo de onde vêm os nomes para os novos atributos
values_from, o atributo de onde vêm os valores dos novos atributos

Note também a função separate(), que divide colunas de acordo com caracteres separadores.

siab_familias <- read_csv2(
    "dados/siab/cadastro_numero_familias.csv", 
    skip = 3,
    locale = locale(encoding = "latin1" ) 
    ) %>% 
    pivot_longer(
        cols = -`Município`,
        names_to = "data",
        values_to = "familias"
    ) %>% 
    rename(
        municipio = `Município`
    ) %>% 
    separate(
        col = municipio,
        into = c("cod_municipio", "municipio"),
        sep = " ",
        extra = "merge"  
    ) %>% 
    mutate(
        cod_municipio = as.integer(cod_municipio),
        data = as.integer(data),
        familias = as.integer(familias)
    )


head(siab_familias)

## # A tibble: 6 x 4
##   cod_municipio municipio              data familias
##           <int> <chr>                 <int>    <int>
## 1        110001 Alta Floresta D'Oeste  1998        2
## 2        110001 Alta Floresta D'Oeste  1999     3458
## 3        110001 Alta Floresta D'Oeste  2000     4668
## 4        110001 Alta Floresta D'Oeste  2001     5078
## 5        110001 Alta Floresta D'Oeste  2002     5078
## 6        110001 Alta Floresta D'Oeste  2003     5077

Pivoteando nosso data frame com muitas colunas (cont.)

Para pegar mais informações dos municipios, vamos ler um arquivo baixado do IBGE

municipios <- read_excel("dados/ibge/populacao.xlsx", skip = 2, col_names = TRUE) 


head(municipios, n = 30)

## # A tibble: 30 x 4
##    Cód.    Município             Ano   `Tabela 6579 - População residente ~
##    <chr>   <chr>                 <chr>                                <dbl>
##  1 1100015 Alta Floresta D'Oest~ 2001                                 26919
##  2 <NA>    <NA>                  2002                                 27237
##  3 <NA>    <NA>                  2003                                 27563
##  4 <NA>    <NA>                  2004                                 29001
##  5 <NA>    <NA>                  2005                                 28629
##  6 <NA>    <NA>                  2006                                 29005
##  7 <NA>    <NA>                  2008                                 24577
##  8 <NA>    <NA>                  2009                                 24354
##  9 <NA>    <NA>                  2011                                 24228
## 10 <NA>    <NA>                  2012                                 24069
## # ... with 20 more rows

Ops… células mescladas

Não deixe o ódio tomar você…

A função fill() pode te ajudar

municipios_ajeitado <- municipios %>%  
  rename(
    cod_municipio = 1,
    nome_municipio = 2,
    ano = 3,
    populacao = 4
  ) %>% 
  fill(cod_municipio, .direction = "down") %>% 
  fill(nome_municipio, .direction = "down") %>% 
  mutate(UF = str_extract(nome_municipio,"\\([A-Z]{2}\\)")) %>%
  mutate(UF = str_extract(UF,"[A-Z]{2}")) %>% 
  mutate(
    cod_municipio = as.integer(cod_municipio),
    ano = as.integer(ano)
  )

head(municipios_ajeitado)

## # A tibble: 6 x 5
##   cod_municipio nome_municipio               ano populacao UF   
##           <int> <chr>                      <int>     <dbl> <chr>
## 1       1100015 Alta Floresta D'Oeste (RO)  2001     26919 RO   
## 2       1100015 Alta Floresta D'Oeste (RO)  2002     27237 RO   
## 3       1100015 Alta Floresta D'Oeste (RO)  2003     27563 RO   
## 4       1100015 Alta Floresta D'Oeste (RO)  2004     29001 RO   
## 5       1100015 Alta Floresta D'Oeste (RO)  2005     28629 RO   
## 6       1100015 Alta Floresta D'Oeste (RO)  2006     29005 RO

Combinando tibbles (dplyr)

Para juntar as informações de dois tibbles em um só, podemos fazer isso de três formas

Combinando tibbles (tidyr): funções de join

Eis as funções de join de data frames

Fonte: (RStudio 2019b)

Combinando tibbles (tidyr): exemplo funções de join:

Anteriormente, pegamos informações do cadastro de famílias…

head(siab_familias)

## # A tibble: 6 x 4
##   cod_municipio municipio              data familias
##           <int> <chr>                 <int>    <int>
## 1        110001 Alta Floresta D'Oeste  1998        2
## 2        110001 Alta Floresta D'Oeste  1999     3458
## 3        110001 Alta Floresta D'Oeste  2000     4668
## 4        110001 Alta Floresta D'Oeste  2001     5078
## 5        110001 Alta Floresta D'Oeste  2002     5078
## 6        110001 Alta Floresta D'Oeste  2003     5077

E de municípios

head(municipios_ajeitado)

## # A tibble: 6 x 5
##   cod_municipio nome_municipio               ano populacao UF   
##           <int> <chr>                      <int>     <dbl> <chr>
## 1       1100015 Alta Floresta D'Oeste (RO)  2001     26919 RO   
## 2       1100015 Alta Floresta D'Oeste (RO)  2002     27237 RO   
## 3       1100015 Alta Floresta D'Oeste (RO)  2003     27563 RO   
## 4       1100015 Alta Floresta D'Oeste (RO)  2004     29001 RO   
## 5       1100015 Alta Floresta D'Oeste (RO)  2005     28629 RO   
## 6       1100015 Alta Floresta D'Oeste (RO)  2006     29005 RO

Os códigos são diferentes, mas

de_para_codigos <- read_csv2(
  "http://blog.mds.gov.br/redesuas/wp-content/uploads/2018/06/Lista_Munic%C3%ADpios_com_IBGE_Brasil_Versao_CSV.csv",
  locale = locale(encoding = "latin1")
  ) %>% 
  select(IBGE, IBGE7)


head(de_para_codigos)

## # A tibble: 6 x 2
##     IBGE   IBGE7
##    <dbl>   <dbl>
## 1 110001 1100015
## 2 110002 1100023
## 3 110003 1100031
## 4 110004 1100049
## 5 110005 1100056
## 6 110006 1100064

Agora juntamos as informações do Siab com as informações dos municípios

siab_familias %>% 
  inner_join(de_para_codigos, by = c("cod_municipio" = "IBGE")) %>%
  inner_join(municipios_ajeitado, by = c("IBGE7" = "cod_municipio", "data" = "ano") ) %>%
  View()

head(siab_familias)

## # A tibble: 6 x 4
##   cod_municipio municipio              data familias
##           <int> <chr>                 <int>    <int>
## 1        110001 Alta Floresta D'Oeste  1998        2
## 2        110001 Alta Floresta D'Oeste  1999     3458
## 3        110001 Alta Floresta D'Oeste  2000     4668
## 4        110001 Alta Floresta D'Oeste  2001     5078
## 5        110001 Alta Floresta D'Oeste  2002     5078
## 6        110001 Alta Floresta D'Oeste  2003     5077

Seria legal saber qual o tamanho médio das famílias

numero_medio_familias <- read_excel("dados/ibge/pessoas_por_familia.xlsx", skip = 2) %>% 
  select(c(1, 3, 7)) %>% 
  rename(
    cod_municipio = 1,
    pessoas = 2,
    numero = 3
  ) %>% 
  
  View()

VISUALIZAÇÃO DE DADOS

NÃO!!! Distorções

Nos próximos slides vemos antipadrões de visualização e dados.

Fonte: (Exame 2018)

NÃO!!! Distorções II

Às vezes mesmo sem querer (será?)

PS. gráfico de 2014

Fonte: (Alves 2014)

NÃO!!! Estilos que atrapalham

Gráfico com sombra, 3D, estilos que dificultam a interpretação

Fonte: (Healy 2018)

NÃO!!! 3D indiscriminado

O pessoal que usa Excel muitas vezes ama 3D, mas…

(Healy 2018)

NÃO!!! Pizza maior que 100%

Fonte: (White 2015)

MELHOR AINDA: Pizza só meio a meio

Assim como os restaurantes devemos servir pizzas de dois sabores no máximo

Tentem descobrir qual a menor e a menor categoria em cada caso

Fonte: (Viz 2018)

NÃO!!! Pizza só meio a meio (cont.)

Fonte: (Viz 2018)

Com muita parcimônia: gráficos de linhas com dois eixos

Gráficos de dois eixos podem mostrar relações espúrias muito facilmente. E elas dependem da escolha da escala e dos limites dos eixos.

Fonte (Rost 2018)

O mesmo dado pode levar aos seguintes gráficos (e suas soluções)

Dicas para uma boa visualização

Algumas dicas para boa visualização de dados:

Enfatize o dado relevante
Integre texto e gráfico
Use fontes e cores diferentes do padrão
Principalmente em colunas e barras, faça o eixo começar de ZERO

loadfonts(device = "win")

by_country <- organdata %>% group_by(consent_law, country) %>%
    summarize_if(is.numeric, list(mean = mean, sd = sd), na.rm = TRUE) %>%
    ungroup()


p <- ggplot(data = by_country,
            mapping = aes(x = gdp_mean, y = health_mean))

p + geom_point( color = "coral4") +
    geom_text_repel(data = subset(by_country, gdp_mean > 25000),
                    mapping = aes(label = country), 
                    size = 3,
                    color = "coral4",
                    family="Bookman Old Style"

                    
                    
                    ) +

  labs(y = "Gastos com saúde per cap.", x = "PIB per cap." ) +
  scale_x_continuous(
    labels = number_format(decimal.mark = ",", big.mark = "."),
    limits = c(0,NA)
  ) +
  scale_y_continuous(
    labels = number_format(decimal.mark = ",", big.mark = "."),
    limits = c(0,NA)
  ) +
  theme_minimal(
  ) +
  ggtitle(
    label = "Gastos em saúde x PIB"
  ) +
  theme(
    text=element_text(family="Bookman Old Style", color = "coral4"),
    axis.text =  element_text(colour = "coral4"),
    panel.background = element_rect(fill = "beige"),
    panel.grid.minor =   element_line(colour = "bisque1"),
    panel.grid.major =  element_line(colour = "bisque3")
  )

Dicas para uma boa visualização

Use small multiples: divida o gráfico em gráficos iguais menores cada um com parte dos dados, seguindo uma regra categórica

#("azure2", "chocolate", "steelblue4", "darkslategray", "slategray3", "slategray4")

p <- ggplot(data = gapminder, mapping = aes(x = year, y = gdpPercap))
p + geom_line(color="chocolate", aes(group = country)) +
    geom_smooth(size = 1.1, method = "loess", se = FALSE, color = "tan4") +
    scale_y_log10(labels=scales::dollar) +
    facet_wrap(~ continent, ncol = 5) +
    labs(x = "Year",
         y = "GDP per capita",
         title = "GDP per capita on Five Continents")+
  theme_minimal() +
  theme(
    text=element_text(family="CMU Serif", color = "darkslategray", size = 8),
    axis.text =  element_text(colour = "darkslategray"),
    panel.background = element_rect(fill = "azure2"),
    panel.grid.minor =   element_line(colour = "slategray2"),
    panel.grid.major =  element_line(colour = "slategray2"),
    strip.text = element_text(family="CMU Serif", colour = "darkslategray"),
    aspect.ratio = 1
  )

GGPLOT2

ggplot2 é a biblioteca de visualização de dados moderna mais usada no R.

Muitos dos problemas listados anteriormente são tratados por ela. É até difícil causar alguns dos problemas anteriores, por exemplo as distorções. É preciso muito malabarismo para produzir um gráfico com 2 eixos.

Por outro lado, a biblioteca facilita muito a criação de small multiples, a criação de gráficos personalizados e complexos

Filosofia

GGPLOT é baseada na filosofia de Tufte (Tufte 1973) e Wilkinson (Wilkinson 1999), em que os gráficos são construídos a partir de dois princípios:

Um gráfico é uma construção em camadas
Os dados ganham sentido a partir de mapeamentos a escalas, que são mapeadas a geometrias

Camadas da GGPLOT2

A GGPLOT2 parece mais complicada de usar do que funções como a plot() da biblioteca base do R.

Talvez porque as pessoas não são sempre apresentadas às camadas da gramática “Grammar of Graphics” que fundamenta a biblioteca.

As camadas da ggplot

Fonte: (Scavetta 2018)

Descrição das camadas

As camadas da ggplot são as seguintes:

Dados (data)
Aesthetics (aes()): mapeamento das colunas dos dados a escalas
Geometries geoms_(): elementos visuais que mostram os dados nas escalas
Facets: small multiples
Statistics: representações do dados transformadas por operações matemáticas
Coordinates
Themes: “non-data ink”, ou seja, os elementos que não são diretamente plotados pela existência de dados

Elementos das camadas

Fonte: (Scavetta 2018)

Montando o gráfico

Usando as camadas, vamos montando o gráfico

Vamos usar para este exemplo um dataset clássico, de espécimes de flores: iris

glimpse(iris)

## Observations: 150
## Variables: 5
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9,...
## $ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1,...
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5,...
## $ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1,...
## $ Species      <fct> setosa, setosa, setosa, setosa, setosa, setosa, s...

Montando o gráfico: data, aes, geom

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +     
  geom_point(alpha = 0.6)

Veja que, pela variável ser representada de forma discreta (uma casa decimal), há sobreposição de pontos. para dar a real impressão de quantos pontos existem, o ideal é inserir um ruído com geom_jitter().

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +     
  geom_jitter(alpha = 0.6)

Montando o gráfico: data, aes, geom, facet

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +     
  geom_jitter(alpha = 0.6) +
  facet_grid(. ~ Species)

Montando o gráfico: data, aes, geom, facet, statistics

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +     
  geom_jitter(alpha = 0.6) +     
  facet_grid(. ~ Species) +
  stat_smooth(method = "lm", se = F, col = "red")

Montando o gráfico: data, aes, geom, facet, statistics, coord

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +     
  geom_jitter(alpha = 0.6) +     
  facet_grid(. ~ Species) +
  stat_smooth(method = "lm", se = F, col = "red") +
  coord_flip()

Montando o gráfico: data, aes, geom, facet, statistics, coord, theme

wes <- wes_palette("Royal2", 5, "discrete")

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +     
  geom_jitter(alpha = 0.6, color = wes[1]) +     
  facet_grid(. ~ Species) +
  stat_smooth(method = "lm", se = F, col = wes[1]) +
  coord_flip() +
  theme_minimal() +
  theme(
    text=element_text(family="Bradley Hand ITC", color = wes[1], size = 16),
    axis.text =  element_text(colour = wes[1]),
    panel.background = element_rect(fill = "white"),
    panel.grid.minor =   element_line(colour = wes[2]),
    panel.grid.major =  element_line(colour = wes[3]),
    strip.text = element_text(family="Bradley Hand ITC", color = wes[1]),
    panel.border = element_rect(colour = wes[4], fill = NA)
  )

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font
## family not found in Windows font database

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x,
## x$y, : font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x,
## x$y, : font family not found in Windows font database

Escolhendo o melhor gráfico

Escolhendo o melhor gráfico: distribuição (variável discreta)

gols <- read_csv("dados/football_events/ginf.csv") %>% 
  select(fthg, ftag) %>% 
  pivot_longer(cols = everything(),  names_to = "casa_fora", values_to = "gols")


ggplot(gols) + 
  geom_histogram(
    aes(x = gols),
    binwidth = 1
  ) +
  scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
  theme_minimal()

result <- fitdistr(  gols$gols  , densfun = "Poisson" ) 

pois <- rpois(length(gols$gols ), lambda = result$estimate ) %>% 
  enframe(value = "gols") %>% 
  mutate(tipo = "poisson")


real_simulado <- gols %>% 
  mutate(tipo = "real" ) %>% 
  bind_rows(pois)

ggplot() + 
  geom_histogram(
    data = real_simulado,
    aes(x = gols, group = tipo, fill = tipo),
    binwidth = 1,
    show.legend = TRUE,
    position = "identity",
    alpha = 0.5
  ) +
  scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
  geom_vline(aes(xintercept = result$estimate)) +
  geom_text(aes(x = result$estimate, y = 7000 ),nudge_x = 0.3,   label = number(result$estimate, accuracy = 0.01 )) +
  theme_minimal()

Escolhendo o melhor gráfico: acumulada empírica

O gráfico mostrando gráfco com a função de probablidade acumulada empírica mostra propriedades que o histograma não mostra

ggplot(gols, aes(x = gols)) +
  stat_ecdf(geom = "step") +
  scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
  theme_minimal()

Escolhendo o melhor gráfico: distribuição de variável contínua

Para variáveis contínuas, faz sentido usar o gráfico de densidade de probabilidade

ufc_pesos <- read_csv("dados/ufc/data.csv") %>% 
  filter(weight_class == "Heavyweight") %>% 
  select(
    R_Height_cms,
    B_Height_cms,
    Winner
  ) %>% 
  transmute(
    altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms  ),
    altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms  )
  ) %>% 
  pivot_longer(cols = everything(), names_to = "resultado", values_to = "altura")


ggplot(ufc_pesos) + 
  geom_density(aes(x = altura)) +
  theme_minimal()

Escolhendo o melhor gráfico: distribuição de variável contínua. ECF

Mais uma vez a distribuição empírica nos deixa ver mais coisa: os quantis

ggplot(ufc_pesos) + 
  stat_ecdf(aes(x = altura), geom = "density" ) +
  theme_minimal()

## Warning: Removed 1 rows containing non-finite values (stat_ecdf).

Escolhendo o melhor gráfico: comparando distribuições

ggplot(ufc_pesos) + 
  geom_density(aes(x = altura, color = resultado)) +
  theme_minimal()

## Warning: Removed 1 rows containing non-finite values (stat_density).

ggplot(ufc_pesos) + 
  stat_ecdf(aes(x = altura, color = resultado), geom = "density" ) +
  theme_minimal()

## Warning: Removed 1 rows containing non-finite values (stat_ecdf).

Escolhendo o melhor gráfico: comparando muitas distribuições

O tipo de gráfico gerado por geom_desity_ridges() ajuda a comparar várias distribuições.

ufc_pesos <- read_csv("dados/ufc/data.csv") %>% 
  select(
    weight_class,
    R_Height_cms,
    B_Height_cms,
    Winner
  ) %>% 
  transmute(
    weight_class,
    altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms  ),
    altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms  )
  ) %>% 
  pivot_longer(cols = -weight_class, names_to = "resultado", values_to = "altura")


ggplot(ufc_pesos) +
  geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
  theme_ridges() +
  scale_x_continuous(breaks = seq(150,by = 5, to = 220))

Legal a ideia, mas não queremos comparar mulheres também

Introduzindo stringr

A biblioteca stringr facilita muito lidar com strings. Por exemplo detectar um pedaço de string em outra.

ufc_pesos <- read_csv("dados/ufc/data.csv") %>% 
  filter(!str_detect(weight_class, "Women")) %>% 
  select(
    weight_class,
    R_Height_cms,
    B_Height_cms,
    Winner
  ) %>% 
  transmute(
    weight_class,
    altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms  ),
    altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms  )
  ) %>% 
  pivot_longer(cols = -weight_class, names_to = "resultado", values_to = "altura")


ggplot(ufc_pesos) +
  geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
  theme_ridges() +
  scale_x_continuous(breaks = seq(150,by = 5, to = 220))

Melhorou, mas ainda á estranho. Melhor seria se as classes de peso estivessem em ordem

Introduzindo forcats

No R, as variáveis categóricas são chamadas de fator. Cada valor possível de um fator é representadas internamente com um identificador numérico, e não com a própria string.

Cada valor possível da variável categórica também é chamada de level

Podemos manipular essas representações facilmente com a forcats.

No caso, queremos ordenar nosso fator com a ordem as classes de peso no UFC

A função fct_relevel possibilita a ordenação manual de uma variável factor

ufc_pesos_ordem <- ufc_pesos %>% 
  mutate(
    weight_class = fct_relevel(weight_class,
      c(
      "Flyweight",
      "Bantamweight",
      "Featherweight",
      "Lightweight",
      "Welterweight",
      "Middleweight",
      "Light Heavyweight",
      "Heavyweight",
      "Catch Weight",
      "Open Weight"
      )
    )  
  )


ggplot(ufc_pesos_ordem) +
  geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
  theme_ridges() +
  scale_x_continuous(breaks = seq(150,by = 5, to = 220))

Relação entre variáveis: scatter plot

Os gráficos do tipo scatter plot possibilitam a visualização de relações entre as variáveis

Vamos brincar um pouco com os dados do Banco Mundial

taxa_homicidios <- wb("VC.IHR.PSRC.P5", country = "all", mrv = 10, POSIXct = TRUE ) %>% 
  select(iso3c, value, date) %>% 
  group_by(iso3c) %>% 
  top_n(1, date) %>% 
  rename(homicidios = value)

gini <- wb("SI.POV.GINI", country = "all", mrv = 10, POSIXct = TRUE ) %>% 
  select(iso3c, value, date) %>% 
  group_by(iso3c) %>% 
  top_n(1, date) %>% 
  rename(desigualdade = value)

pib_per_capita <- wb("NY.GDP.PCAP.PP.KD", country = "all", mrv = 10, POSIXct = TRUE ) %>% 
  select(iso3c, value, date) %>% 
  group_by(iso3c) %>% 
  top_n(1, date) %>% 
  rename(riqueza = value)

pib_gini_homi <- taxa_homicidios %>% 
  inner_join(gini, by = c("iso3c")) %>% 
  inner_join(pib_per_capita, by = c("iso3c")) %>% 
  left_join(codelist, by = c("iso3c") ) %>% 
  select(riqueza, desigualdade, homicidios, continent, region) %>% 
  pivot_longer(cols = c("riqueza", "desigualdade"), names_to = "tipo_indice", values_to = "indice" ) %>% 
  mutate(tipo_indice = str_to_title(tipo_indice))

ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
  geom_point(alpha = 0.4) +
  facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
  geom_smooth() +
  theme_minimal() +
  labs(x = "", y = "Homicídios por 100 mil") +
  scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = "."))

Relação entre variáveis: scatter plot (cont.)

É possível visualizar os eixos em log na base 10

ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
  geom_point(alpha = 0.4) +
  facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
  geom_smooth() +
  theme_minimal() +
  labs(x = "", y = "Homicídios por 100 mil") +
  scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
  scale_y_log10()

Relação entre variáveis: scatter plot (cont.)

É possível também usar as cores para olhar uma terceira característica

ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
  geom_point(alpha = 0.4, aes(color = continent )) +
  facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
  geom_smooth() +
  theme_minimal() +
  labs(x = "", y = "Homicídios por 100 mil") +
  scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
  scale_y_log10()  +
  theme(legend.position = "top")

A linha de regressão também pode ser estimada para cada grupo. Neste caso vamos usar o método de regressão linear, pois o LOESS original nao funciona nada bem com muito poucos dados (nem a linear, mas…)

ggplot(pib_gini_homi, aes(y = homicidios, x = indice, color = continent)) +
  geom_point(alpha = 0.4 ) +
  facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
  geom_smooth(se = FALSE, method = "lm") +
  theme_minimal() +
  labs(x = "", y = "Homicídios por 100 mil") +
  scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
  scale_y_log10()

analise_desigualdade <- pib_gini_homi %>% 
  filter(tipo_indice == "Desigualdade")

ggplot(analise_desigualdade, aes(y = homicidios, x = indice, color = continent)) +
  geom_point(alpha = 0.4) +
  facet_wrap(~tipo_indice, ncol = 1) +
  geom_smooth(se = FALSE, method = "lm") +
  theme_minimal() +
  labs(x = "Índice de Gini", y = "Homicídios por 100 mil") +
  scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
  scale_y_log10() +
  facet_wrap(~continent)  +
  theme(legend.position = "top")

ggplot(analise_desigualdade, aes(y = homicidios, x = indice, color = continent)) +
  geom_text(aes(label = iso3c), size = 2) +
  facet_wrap(~tipo_indice, ncol = 1) +
  geom_smooth(se = FALSE, method = "lm") +
  theme_minimal() +
  labs(x = "Índice de Gini", y = "Homicídios por 100 mil") +
  scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
  scale_y_log10() +
  facet_wrap(~continent) +
  theme(legend.position = "top")

Relação entre variáveis: Três variáveis contínuas

Podemos avaliar três variáveis contínuas usando um gráfico

Primeiro vamos colocar em wide as duas variáveis que tínhamos deixado long.

pib_gini_homi_wide <- pib_gini_homi %>% 
  pivot_wider(names_from = tipo_indice, values_from = indice )
  

head(pib_gini_homi_wide)

## # A tibble: 6 x 6
## # Groups:   iso3c [6]
##   iso3c homicidios continent region                    Riqueza Desigualdade
##   <chr>      <dbl> <chr>     <chr>                       <dbl>        <dbl>
## 1 ALB          2.3 Europe    Southern Europe            12316.         29  
## 2 DZA          1.4 Africa    Northern Africa            13737.         27.6
## 3 AGO          4.8 Africa    Middle Africa               5725.         42.7
## 4 ARG          5.1 Americas  South America              18288.         41.2
## 5 ARM          2.4 Asia      Western Asia                9178.         33.6
## 6 AUS          0.8 Oceania   Australia and New Zealand  45378.         35.8

A escala de cor usada abaixo, Viridis, oferece a mesma sensação para colorido e preto e branco e é feita para ajudar daltônicos.

O log foi usado aqui para evitar que os valores extremos dificultem a visualização da escala de cores.

ggplot(pib_gini_homi_wide) +
  geom_point(aes(color = homicidios, x = Desigualdade , y = Riqueza )) +
  scale_color_viridis_c(trans = "log10", direction = -1, option = "inferno") +
  theme_minimal()

Relação entre variáveis: Três variáveis contínuas

hdi <- read_csv("dados/hdi/atlas.csv") %>% 
  rename(municipio = `município` ) %>% 
  select(
    ano,
    rdpc,
    gini,
    municipio,
    uf
  ) %>% 
  filter(ano == 2010)
  

violencia <- extract_text("dados/atlasviolencia/8099-tabelamunicipiostodossite.pdf", encoding = "UTF-8") %>% 
  str_split(pattern = fixed("\n")) %>% 
  unlist() %>% 
  enframe( value = "linha") %>% 
  mutate(
    UF = str_extract(linha, "[A-Z]{2}") ,
    municipio = str_extract(linha, "(?<=[A-Z]{2} ).+?(?=[0-9])"),
    numeros = str_extract_all(linha, "(?<= )[0-9 ,.]*") 
  ) %>% 
  unnest_legacy() %>% 
  filter(str_length(numeros)>1) %>% 
  mutate(
    numeros = str_trim(numeros),
    municipio = str_trim(municipio) %>% str_to_upper()
  ) %>% 
  separate(
    col = numeros,
    into = c("populacao", "homicidios", "ocultos", "taxa_homicidios"),
    sep = " "
  ) %>% 
  mutate(taxa_homicidios = parse_number(taxa_homicidios, locale = locale(decimal_mark = ",")) )
  

ufs <- municipios_ajeitado %>% 
  select(
    UF,
    cod_municipio
  ) %>% 
  mutate(
    cod_uf = cod_municipio %/% 100000
  ) %>% 
  select(-cod_municipio) %>% 
  distinct() %>% 
  filter(!is.na(cod_uf)) 
  
hdi %<>% inner_join(ufs, by = c("uf" = "cod_uf"), suffix = c("_old", "")) 


hdi_violencia <- hdi %>% 
  inner_join(violencia, by = c("municipio", "UF") ) %>% 
  mutate(
    regiao = case_when(
      uf %/% 10 == 1 ~ "Norte",
      uf %/% 10 == 2 ~ "Nordeste",
      uf %/% 10 == 3 ~ "Sudeste",
      uf %/% 10 == 4 ~ "Sul",
      uf %/% 10 == 5 ~ "Centro-Oeste",
      TRUE ~ NA_character_
      
    ) 
  )

ggplot(hdi_violencia) +
  geom_jitter(aes(color = taxa_homicidios, x = gini , y = rdpc ), alpha = 0.2   ) +
  scale_color_gradient(low = "blue", high = "red", trans = "log10" ) +
  facet_wrap(~regiao) + 
  theme_minimal()

## Warning: Transformation introduced infinite values in discrete y-axis

Relação entre variáveis: Duas discretas e uma contínua

ufc_data <- read_csv("dados/ufc/data.csv") %>% 
  select(weight_class, no_of_rounds )

ufc_raw_data <- read_csv2("dados/ufc/raw_total_fight_data.csv") %>% 
  select(last_round)


ufc_tudo <- bind_cols(ufc_data, ufc_raw_data) %>% 
  filter(no_of_rounds == 3) %>% 
  group_by(weight_class) %>% 
  mutate(n_classe = n()) %>%
  group_by(weight_class, last_round) %>% 
  summarise(n_round = n(), n_classe = mean(n_classe)) %>% 
  mutate(frac_lutas = n_round/n_classe ) %>% 
  ungroup() %>% 
  filter( weight_class %in%
      c(
      "Flyweight",
      "Bantamweight",
      "Featherweight",
      "Lightweight",
      "Welterweight",
      "Middleweight",
      "Light Heavyweight",
      "Heavyweight"
            
          )) %>% 
  mutate(
    weight_class = fct_relevel(weight_class,
      c(
      "Flyweight",
      "Bantamweight",
      "Featherweight",
      "Lightweight",
      "Welterweight",
      "Middleweight",
      "Light Heavyweight",
      "Heavyweight"
      )
    )
  )

O gráfico gerado por geom_tile() ajuda a visualizar duas variáveis discretas ordinais e uma contínua.

ggplot(ufc_tudo, aes(x = last_round, y = weight_class)) +
  geom_tile(aes(fill = frac_lutas )) +
  geom_text(aes(label = percent(frac_lutas))   ) +
  scale_fill_gradient(low = "white", high = "darkgreen", labels = percent ) +
  theme_minimal()

Relação entre variáveis: Várias variáveis

hdi_desc <-  read_csv("dados/hdi/desc.csv")

hdi_tudo <-  read_csv("dados/hdi/atlas.csv") %>% 
  mutate(
    regiao = case_when(
      uf %/% 10 == 1 ~ "N",
      uf %/% 10 == 2 ~ "NE",
      uf %/% 10 == 3 ~ "SE",
      uf %/% 10 == 4 ~ "S",
      uf %/% 10 == 5 ~ "CO",
      TRUE ~ NA_character_
    )
  ) %>% 
  select(
      espvida,
      anos_estudo = e_anosestudo,
      gini,
      pind,
      pren10ricos,
      prentrab,
      t_banagua,
      regiao
  )


ggpairs(hdi_tudo)  +
  theme_minimal()

Composição: no tempo, poucos períodos, valores relativos

A função geom_col() é melhor para poucos períodos

Pelo mapeamentos no eixo x e y, haveria valores no mesmo lugar do eixo cartesiano.

Para sanar esse conflito, usamos o position - parâmetro da geom_col().

Pata o caso em que queremos comparar valores relativamente a cada período, devemos empilhar as categorias de forma que a coluna tenha o mesmo tamanho a cada período.

Para isso, usamos o valor fill para o parâmetro position

jogos <- read_csv("dados/football_events/ginf.csv") %>% 
  mutate(
    resultado = case_when(
      fthg > ftag ~ "Casa",
      fthg < ftag ~ "Fora",
      TRUE ~ "Empate"
      )
  ) %>% 
  count(season, country, resultado) 



ggplot(jogos) + 
  geom_col(aes(y = n, fill = resultado, x = as.factor(season)), position = "fill") +
  facet_wrap(~country) +
  theme_minimal() +
  scale_fill_manual(values = wes_palette("Royal2")) +
  scale_y_continuous(labels = percent) +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "top",
    text = element_text(family = "Rockwell Condensed")
  ) +
  labs(x = "Temporada", y = "% Vitórias", fill = "Resultado")

Composição: no tempo, poucos períodos, valores absolutos

Para compararmos de forma que o total de cada período apareça, em valores absolutos, devemos usar o valor stack

jogos <- read_csv("dados/football_events/ginf.csv") %>% 
  filter(season != 2017) %>% 
  select(season, country, Fora = ftag, Casa = fthg) %>% 
  pivot_longer(cols = c(Fora, Casa), names_to = "Time", values_to = "Gols"  ) %>% 
  group_by(season, country, Time) %>% 
  summarise(Gols = sum(Gols)) 


ggplot(jogos) + 
  geom_col(aes(y = Gols, fill = Time, x = as.factor(season)), position = "stack") +
  facet_wrap(~country) +
  theme_minimal() +
  scale_fill_manual(values = wes_palette("Royal2")) +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "top",
    text = element_text(family = "Rockwell Condensed")
  ) +
  labs(x = "Temporada", y = "Gols", fill = "Time")

Composição: no tempo, muitos períodos, valores relativos

Para muito períodos o ideal é usar áreas empilhadas, com a função geom_area()

Repare que há muitos países, ou seja, muitas categorias.

Usando a biblioteca forcats e sua função fct_lump() conseguimos criar uma categoria de “Outros”

gdps <- wb(indicator = "NY.GDP.MKTP.PP.KD",  country = "countries_only") 


gdps_tratado <- gdps %>% 
  mutate(date = as.integer(date) ) %>% 
  group_by(country) %>% 
  mutate(ultimo_gdp = last(value, order_by = date)) %>% 
  ungroup() %>%  
  mutate(country = fct_lump(country, n = 7, w = ultimo_gdp, other_level = "Outros")) %>% 
  group_by(country, date) %>% 
  summarise(value = sum(value),
            ultimo_gdp = sum(ultimo_gdp)) %>% 
  ungroup() %>% 
  mutate(country = fct_reorder(country, ultimo_gdp ))


ggplot(gdps_tratado ) +
  geom_area(aes(x = date, y = value, fill = country), position = "fill") +
  theme_minimal() +
  scale_fill_brewer(palette = "Accent") +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "top",
    text = element_text(family = "CMU Classical Serif")
  ) +
  labs(x = "Ano", y = "PIB PPP", fill = "País")

Composição: no tempo, muitos períodos, valores absolutos

gdps_tratado <- gdps %>% 
  mutate(date = as.integer(date) ) %>% 
  group_by(country) %>% 
  mutate(ultimo_gdp = last(value, order_by = date)) %>% 
  ungroup() %>%  
  mutate(country = fct_lump(country, n = 7, w = ultimo_gdp, other_level = "Outros")) %>% 
  group_by(country, date) %>% 
  summarise(value = sum(value),
            ultimo_gdp = sum(ultimo_gdp)) %>% 
  ungroup() %>% 
  mutate(country = fct_reorder(country, ultimo_gdp ))


ggplot(gdps_tratado ) +
  geom_area(aes(x = date, y = value, fill = country) ) +
  theme_minimal() +
  scale_fill_brewer(palette = "Accent") +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "top",
    text = element_text(family = "CMU Serif Extra")
  ) +
  labs(x = "Ano", y = "PIB PPP", fill = "País")

Composição: estática, valores relativos

Aqui cabe uma torta se houver poucas categorias

Para gerar um gráfico de torta é preciso gerar um gráfico de barras e tornar a coordenada polar com uso de coord_polar

gdps_tratado <- gdps %>% 
  mutate(date = as.integer(date)) %>% 
  filter(date == max(date) ) %>% 
  mutate(country = fct_lump(country, n = 1, w =  value, other_level = "Resto")) %>% 
  group_by(country) %>% 
  summarise(PIB = sum(value))

  
ggplot(gdps_tratado, aes(x = "", y = PIB, fill = country)) +
  geom_bar( width = 1, stat = "identity", color = "white") +
  coord_polar("y", start = 0) +
  scale_fill_brewer(palette = "Accent") +
  theme_void()

Composição: estática, valores relativos, muitos

gdps_tratado <- gdps %>% 
  mutate(date = as.integer(date)) %>% 
  filter(date == max(date) ) %>% 
  mutate(country = fct_lump(country, n = 50, w =  value, other_level = "Resto")) %>% 
  group_by(country) %>% 
  summarise(PIB = sum(value))

  
ggplot(gdps_tratado, aes( fill = country, area = PIB, label = country)) +
  geom_treemap() +
  geom_treemap_text(grow = TRUE) +
  theme(
    legend.position = "none"
  )

Visualização Geoespacial

A biblioteca sf tem várias funções para lidar com dados geoespaciais.

A função st_read lê vários formatos de dados geoespaciais para um objeto da classe sf, ou simple features.

gd <- read_sf("dados/gd/Geração_Distribuída.shp")

Classe sf

sf é uma classe que extende a classe data.frame contendo uma coluna de dados geoespaciais.

Cada linha do data.frame é relacionada a um dado geoespacial, que pode ser um ponto, um polígono etc.

str(gd)

## Classes 'sf', 'tbl_df', 'tbl' and 'data.frame':  180674 obs. of  20 variables:
##  $ X         : num  -39.9 -40.6 -37.1 -48.4 -48.4 ...
##  $ Y         : num  -18.72 -19.38 -11.04 -1.36 -1.38 ...
##  $ USER_IdeAg: num  380 381 6587 371 371 ...
##  $ USER_SigAg: chr  "EDP ES" "ELFSM" "ESE" "CELPA" ...
##  $ USER_CodUn: chr  "196457" "33243" "703489" "7744250" ...
##  $ USER_CodGD: chr  "GD.ES.000.034.069" "GD.ES.000.094.399" "GD.SE.000.093.608" "GD.PA.000.108.959" ...
##  $ USER_DscCl: chr  "Residencial" "Rural" "Residencial" "Residencial" ...
##  $ USER_DscGr: chr  "B1" "B2" "B1" "B1" ...
##  $ USER_DscMo: chr  "Autoconsumo remoto" "Autoconsumo remoto" "Geracao na propria UC" "Geracao na propria UC" ...
##  $ USER_QtdUC: num  2 2 1 1 1 1 1 1 1 1 ...
##  $ USER_CodMu: num  3204906 3203353 2800308 1500800 1500800 ...
##  $ USER_NomMu: chr  "São Mateus" "Marilândia" "Aracaju" "Ananindeua" ...
##  $ USER_SigUF: chr  "ES" "ES" "SE" "PA" ...
##  $ USER_NomRe: chr  "Sudeste" "Sudeste" "Nordeste" "Norte" ...
##  $ USER_SigTi: chr  "UFV" "UFV" "UFV" "UFV" ...
##  $ USER_DscCo: chr  "Radiação solar" "Radiação solar" "Radiação solar" "Radiação solar" ...
##  $ USER_MdaPo: num  4.6 3 3.96 3 2.4 ...
##  $ USER_NomAg: chr  "ESPÍRITO SANTO DISTRIBUIÇÃO DE ENERGIA S.A." "EMPRESA LUZ E FORÇA SANTA MARIA S/A" "ENERGISA SERGIPE - DISTRIBUIDORA DE ENERGIA S.A" "CENTRAIS ELÉTRICAS DO PARÁ S.A. - CELPA" ...
##  $ DthConexao: Date, format: "2018-07-24" "2019-05-21" ...
##  $ geometry  :sfc_POINT of length 180674; first list element:  'XY' num  -39.9 -18.7
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "X" "Y" "USER_IdeAg" "USER_SigAg" ...

Operações espaciais com sf

enderecos <- tribble(
  
  ~nome,    ~endereco,
  "RB1",    "Avenida Rio Branco 1, Centro, Rio de Janeiro",
  "Marques dos Reis",    "Praça Pio X 54, Centro, Rio de Janeiro"
  
) 


enderecos_com_geo <- enderecos %>% 
  mutate_geocode(location = endereco)

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Avenida+Rio+Branco+1,+Centro,+Rio+de+Janeiro&key=xxx

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Pra%C3%A7a+Pio+X+54,+Centro,+Rio+de+Janeiro&key=xxx

enderecos_id <- enderecos_com_geo %>% 
  mutate(id = row_number()) 

comb_1 <- enderecos_id %>% 
  rename_all(function(x){str_glue("{x}_1")})

comb_2 <- enderecos_id %>% 
  rename_all(function(x){str_glue("{x}_2")})



ponto_4326 <- function (lat, lon){
  st_point(c(lat, lon)) %>% 
    st_sfc(crs = 4326)
}

enderecos_combinados <- comb_1 %>% 
  crossing(comb_2) %>% 
  filter(id_1 < id_2) %>% 
  mutate(
    ponto_1 = map2(.x = lon_1, .y = lat_1, .f = ponto_4326 ) ,
    ponto_2 = map2(.x = lon_2, .y = lat_2, .f = ponto_4326 ) ,
    dist = map2(.x = ponto_1, .y = ponto_2, st_distance ) 
  )

Biblioteca ggmap

A biblioteca ggmap deixa que usemos um background em um plot de mapa.

Ele funciona em camadas como o ggplot, e aceita as mesmas funcionalidades para plotar os dados, usando o eixo x como latitude e o eixo y como longitude

No código abaixo, usamos get_map para pegar o fundo, passando os parâmetros desejados. Existem várias fontes de dados e vários tipos de fundo

O tema incluído por theme_inset é mais adequado a mapas: não mostra os eixos, com ticks e rótulos, por exemplo.

As coordenadas são mantidas fixas, também, com uso de coord_fixed

map <- get_map(location = "Brazil",  zoom = 4, source = "google", maptype = "toner-lite")

ggmap(map) +
  coord_fixed() +
  theme_inset()

Plotando objetos no mapa

Podemos plotar os objetos como pontos no mapa, por exemplo, usando geom_point()

ggmap(map) + 
  geom_point(
    data = gd, 
    alpha = 0.01,
    size = 0.5,
    aes(
      x = X,
      y = Y
    ),
    color = "darkblue"
  ) +
  theme_inset()

Plotando a densidade de objetos no mapa

Usando stat_density_2d, podemos plotar

map <- get_map(location = "Brazil", zoom = 4,  source = "google", maptype = "toner-lite")

## maptype = "toner-lite" is only available with source = "stamen".

## resetting to source = "stamen"...

## Source : https://maps.googleapis.com/maps/api/staticmap?center=Brazil&zoom=4&size=640x640&scale=2&maptype=terrain&key=xxx

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Brazil&key=xxx

## Source : http://tile.stamen.com/toner-lite/4/4/7.png

## Source : http://tile.stamen.com/toner-lite/4/5/7.png

## Source : http://tile.stamen.com/toner-lite/4/6/7.png

## Source : http://tile.stamen.com/toner-lite/4/4/8.png

## Source : http://tile.stamen.com/toner-lite/4/5/8.png

## Source : http://tile.stamen.com/toner-lite/4/6/8.png

## Source : http://tile.stamen.com/toner-lite/4/4/9.png

## Source : http://tile.stamen.com/toner-lite/4/5/9.png

## Source : http://tile.stamen.com/toner-lite/4/6/9.png

ggmap(ggmap = map) +
  stat_density_2d(
    aes(x = X, y = Y, fill = ..level.. ), 
    bins = 70,
    geom = "polygon", 
    data = gd  ,
    alpha = .2
  ) +
  scale_fill_gradientn(colors = brewer.pal(10, "YlOrRd")) +  
  theme_inset()

## Warning in brewer.pal(10, "YlOrRd"): n too large, allowed maximum for palette YlOrRd is 9
## Returning the palette you asked for with that many colors

MODELOS DE APRENDIZADO ESTATÍSTICO

Rudimentos de Aprendizado Estatístico

Aqui apresentaremos RUDIMENTOS de Aprendizado Estatístico com o objetivo de desmistificar alguns conceitos e apresentar algumas considerações que nos ajudarão a iniciar um processo deste tipo.

Desmistificando Aprendizado Estatístico (ou Machine Learning)

Pesquisando imagens no Google vemos as percepções a respeito de “Machine Learning” e “Statistical Learning”

O processo de “Data Science”

(Mello 2019)

Mecânico ou Piloto ?

O processo de Aprendizado Estatístico

Temos acesso a um pedaço dos dados existentes, que usamos para nos ajudar a especificar e treinar o nosso modelo.

Podemos adotar como premissa que existe uma função que traduz como o universo funciona. Esta função determina qual valor \(y\) da variável dependente acontece quando as variáveis independentes assumem valores \(x\), onde \(x\) é um vetor de tamanho \(p\). \(p\) é o número de variáveis dependentes usadas. Esta função não determina perfeitamente \(y\), então temos sempre um \(\epsilon\), um ruído que pode ser devido a variáveis dependentes desconhecidas ou efeitos que não podem ser mensurados.

\[y = f(x) + \epsilon\]

Onde:

\(f(x)\) é uma função desconhecida e \(E(\epsilon) = 0\),

\(\epsilon\) é independente de \(x\) e \(Var(\epsilon) = \sigma^2\) é constante em relação a \(x\).

Nós não temos acesso a esta \(f()\) que representa a VERDADE, mas tentamos estimar uma \(\hat{f}()\)

Fábrica de função que aproxima a realidade

Nosso modelo de Aprendizado Estatístico é uma fábrica de \(\hat{f}()\) que tenta aproximar a \(f()\) verdadeira.

Existem fábricas para todos os gostos. Desde fábricas que produzem funções simples e inteligíveis até funções do tipo caixa preta.

Esta fábrica de \(\hat{f}()\) recebe como insumo variáveis retiradas do ambiente. Essas variáveis podem ser tratadas de forma a deixar as coisas mais fáceis para o modelo, de modo a deixá-lo na cara do gol como Gérson fazia na copa de 70. Esta etapa de preparar as entradas se chama Feature Engineering. Para fazer isso, usamos tudo que aprendemos sobre manipulação de dados.

Modelo treinado

Com o(s) modelo(s) treinado(s), faremos predições a respeito de outros dados futuros ou cuja variável dependente ainda desconhecemos.

Podemos também fazer inferências a respeito de como as variáveis independentes afetam a variáveis dependente.

Inferência x Predição

A predição é o ato de descobrir \(y\) a partir de \(x\).

A inferência é o entendimento de como valores diferentes de x afetam y

Predição

Podemos dividir o erro de previsão em redutível e irredutível.

Temos que:

\[\hat{y} = \hat{f}(x) \]

Então:

\[E(y - \hat{y})^2 = [f(x) - \hat{f}(x)]^2 + Var(\epsilon)\]

Onde

\([f(x) - \hat{f}(x)]^2\) é um erro redutível: podemos sempre estimar uma \(\hat{f}()\) melhor

\(Var(\epsilon)\) é um erro irredutível: não importa quão bem estimemos \(\hat{f}()\) o resíduo \(\epsilon\) está lá

(Hastie, Tibshirani, and Friedman 2001)

Viés e Variância

Podemos dividir o erro redutível em viés e variância.

Se rodarmos \(\hat{f}(x_0)\) treinando o modelo com vários conjuntos de treinamento e testássemos ele com uma entrada específica qualquer para a qual \(y_0 = f(x_0) + \epsilon\):

\[E[y_o - \hat{f}(x_0)] = Var(\hat{f}(x_0)) + [Bias(\hat{f}(x_0))]^2 + Var(\epsilon)\]

\(Var(\hat{f}(x_0)\) é a variância do modelo: o quanto ele muda quando treinamos ele com outro conjunto de treinamento

\([Bias(\hat{f}(x_0))]^2\) é o erro introduzido por tentarmos aproximar um problema complexo da vida real com uso de um modelo simplificado, também chamado de viés do modelo.

Trade-off viés x variância

Os modelos variam muito com relação a “transparência da caixa” e essa relação tem correlação com viés e variância.

Existe um trade-off entre viés e variância, que são medidas que a gente não observa diretamente.

Há uma relação inversamente proporcional (grosso modo) entre:

Flexibilidade do modelo, ou poder de se adaptar a uma relação complexa entre entrada e saída. Modelos com maior flexibilidade tendem a ter viés menor e variância maior.
Interpretabilidade, ou facilidade de o usuário do modelo (piloto) inferir as relações entre as variáveis independentes e as variáveis dependentes. Modelos mais fáceis de interpretar tendem a ser menos poderosos portanto tem viés maior. Mas também variam menos.

(James et al. 2013)

Ponto ótimo no aprendizado e overfitting

Além da falta de interpretabilidade, outra questão pode contar contra os modelos mais complexos, com mais poder de aprender nuances sobre os dados: eles podem aprender características específicas da amostra de treinamento que não podem ser generalizados para o resto dos dados.

Existe um ponto ótimo no poder de aprendizado para que o modelo atinja a melhor performance possível na população como um todo.

Além desse ponto ótimo, onde o erro no conjunto de teste é mínimo, o erro no conjunto de treinamento continua a diminuir, mas o aumento do erro no conjunto de teste mostra que ele perde o poder de generalização. Neste região diz-se que está acontecendo o fenômeno do overfitting.

(Hastie, Tibshirani, and Friedman 2001)

Treino, Validação e Teste

O procedimento mais usado para preparar um modelo para a vida real envolve dividir a nossa base em três pedaços:

Dados de treinamento: que serão usados para treinar o modelo e não servem para avaliar seu desempenho
Dados de validação: servem para comparar modelos de forma que nos possibilite escolher UM MODELO COM UMA ESPECIFICAÇÃO DE HIPERPARÂMETROS (número de nós da rede neural, )
Dados de teste: são separados logo no início do processo, e só são usados numa avaliação final do único modelo escolhido.

Model Selection e Model Assesment

Model Selection: estimação da performance de vários modelos a fim de escolher o melhor

Model Assessment: tendo escolhido o melhor modelo, estimação do erro de generalização em novos dados

K-Fold Cross validation

Uma forma de usar todos os dados (EXCETO OS DADOS DE TESTE) tanto para treinar quando para validar é revezando que partes dos dados estão sendo usados para treinamento e para validação.

Este esquema é chamado de k-fold cross validation

Podemos dividir os dados em \(k\) partes. Cada vez uma parte é escolhida para ser a validação.

(Robot 2019)

Execução de modelos com broom - explorando os dados

Primeiro podemos fazer nossa exploração preliminar.

Como exemplo, podemos executar um pequeno tratamento das features usando a função pivot_longer() de um jeito diferente, só possível na última versão da biblioteca tidyr.

ufc_data <- read_csv("dados/ufc/data.csv") 

ufc_raw_data <- read_csv2("dados/ufc/raw_total_fight_data.csv") 

ufc_raw_fighter <- read_csv2("dados/ufc/raw_fighter_details.csv") 

ufc_raw_pre <- read_csv2("dados/ufc/preprocessed_data.csv") 


ufc_cada_lutador <- bind_cols(ufc_data, ufc_raw_data) %>% 
  pivot_longer(
    cols = matches("[BR]_.*"),
    names_to = c("lutador",".value") ,
    names_sep = "_"
  ) %>% 
  mutate(
    ganhou = str_sub(Winner,1,1) == lutador
  ) %>% 
  mutate(aprov = wins/(wins+losses) )
  
  
ggplot(ufc_cada_lutador) +
  geom_boxplot(aes(x = ganhou, y = aprov )) + 
  facet_wrap(~weight_class) +
  theme_minimal()

ufc_aprend <- bind_cols(ufc_data, ufc_raw_data) %>%
  filter(Winner != "Draw") %>%
  mutate(
    R_aprov = R_wins / (R_wins + R_losses),
    B_aprov = B_wins / (B_wins + B_losses),
  ) %>% 
  select(
    Winner,
    weight_class,
    B_Weight_lbs,
    B_age,
    R_Weight_lbs,
    R_age,
    B_age,
    R_aprov,
    B_aprov

  ) %>%
  mutate(
    winner_color = Winner,
    zebra_blue = if_else(Winner == "Blue",1,0),
    red_mais_velho = R_age - B_age,
    idade_red = R_age
  ) %>% 
  filter(
    !is.na(idade_red) & !is.na(red_mais_velho) & !is.na(zebra_blue)
  )

Estudando algumas features para nosso exemplo

Como exemplo, podemos observar a relação entre a diferença de idade dos lutadores e o vencedor da luta.

Normalmente o lutador vermelho é o favorito. mas podemos ver que quando o lutador vermelho é velho e mais velho que o azul, ele costuma perder.

ggplot(ufc_aprend ) +
  geom_jitter(aes( y = red_mais_velho, x = R_age, color = winner_color ), alpha = 0.15) +
  scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
  facet_wrap( ~weight_class ) +
  theme_minimal()

Estudando algumas features para nossso exemplo

Podemos tentar uma regressão linear para avaliar isso:

\[ zebra = \alpha + \beta_{idade} idade\_vermelho + \beta_{dif} dif\_vermelho\_mais\_velho + \epsilon \]

Podemos perceber que os betas são significativos.

modelo_lm <- lm(zebra_blue ~ idade_red  + red_mais_velho , data = ufc_aprend)

summary(modelo_lm)

## 
## Call:
## lm(formula = zebra_blue ~ idade_red + red_mais_velho, data = ufc_aprend)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6314 -0.3407 -0.2586  0.5898  0.8888 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.051147   0.061605   0.830    0.406    
## idade_red      0.009233   0.002089   4.420 1.01e-05 ***
## red_mais_velho 0.010876   0.001651   6.589 4.91e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4609 on 4865 degrees of freedom
## Multiple R-squared:  0.03404,    Adjusted R-squared:  0.03365 
## F-statistic: 85.73 on 2 and 4865 DF,  p-value: < 2.2e-16

Facilitando a análise de um modelo com `broom`

A biblioteca broom ajuda a extrair informações de vários tipos de modelo com as mesmas funções.

modelo_tidy <- tidy(modelo_lm)

modelo_aumentado <- augment(modelo_lm)

modelo_tidy

## # A tibble: 3 x 5
##   term           estimate std.error statistic  p.value
##   <chr>             <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)     0.0511    0.0616      0.830 4.06e- 1
## 2 idade_red       0.00923   0.00209     4.42  1.01e- 5
## 3 red_mais_velho  0.0109    0.00165     6.59  4.91e-11

head(modelo_aumentado)

## # A tibble: 6 x 10
##   zebra_blue idade_red red_mais_velho .fitted .se.fit .resid    .hat .sigma
##        <dbl>     <dbl>          <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>
## 1          0        32              1   0.357 0.00806 -0.357 3.06e-4  0.461
## 2          0        31             -1   0.326 0.00819 -0.326 3.16e-4  0.461
## 3          0        35             -1   0.363 0.0146  -0.363 1.00e-3  0.461
## 4          1        29              3   0.352 0.00839  0.648 3.31e-4  0.461
## 5          1        26             -6   0.226 0.0103   0.774 5.03e-4  0.461
## 6          0        28             -5   0.255 0.00973 -0.255 4.46e-4  0.461
## # ... with 2 more variables: .cooksd <dbl>, .std.resid <dbl>

Matriz de confusão com a `caret`

A biblioteca caret também ajuda muito.

A função confusionMatrix monta uma matriz de confusão, que mostra estão acontecendo os falsos positivos, falsos negativos, verdadeiros positivos e verdadeiros negativos.

modelo_aumentado_factor <- modelo_aumentado %>% 
  mutate(
    previsao_zebra_blue = as_factor(round(.fitted)),
    zebra_blue = as_factor(zebra_blue)
  ) 

confusionMatrix(modelo_aumentado_factor$previsao_zebra_blue, modelo_aumentado_factor$zebra_blue )

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 3227 1520
##          1   53   68
##                                         
##                Accuracy : 0.6769        
##                  95% CI : (0.6635, 0.69)
##     No Information Rate : 0.6738        
##     P-Value [Acc > NIR] : 0.3293        
##                                         
##                   Kappa : 0.035         
##                                         
##  Mcnemar's Test P-Value : <2e-16        
##                                         
##             Sensitivity : 0.98384       
##             Specificity : 0.04282       
##          Pos Pred Value : 0.67980       
##          Neg Pred Value : 0.56198       
##              Prevalence : 0.67379       
##          Detection Rate : 0.66290       
##    Detection Prevalence : 0.97514       
##       Balanced Accuracy : 0.51333       
##                                         
##        'Positive' Class : 0             
##

Rodando o modelo e plotando um grid de predições usando a biblioteca `purr`

Preparando um grid das variáveis independentes

red_mais_velho <-  seq(
    from = min(ufc_aprend$red_mais_velho, na.rm = TRUE), 
    to =  max(ufc_aprend$red_mais_velho, na.rm = TRUE), 
    length.out = 50
  ) %>% 
  enframe(name = "1", value = "red_mais_velho")

idade_red <-  seq(
    from = min(ufc_aprend$idade_red, na.rm = TRUE), 
    to =  max(ufc_aprend$idade_red, na.rm = TRUE), 
    length.out = 50
  ) %>% 
  enframe(name = "2", value = "idade_red")

weight_classes <- ufc_aprend %>% 
  select(weight_class) %>% 
  distinct()


grid_previsao <- red_mais_velho %>% 
  crossing(idade_red) %>% 
  crossing(weight_classes)


ggplot(grid_previsao %>% filter(weight_class == "Heavyweight")) +
  geom_point(aes(x = idade_red, y = red_mais_velho ), size = 0.01) +
  theme_minimal()

Rodando o modelo e plotando um grid de predições usando a biblioteca `purr`

Rodando o modelo com o grid

previsoes <- augment(modelo_lm, newdata = grid_previsao  ) 


#idade_red  + red_mais_velho

ggplot(ufc_aprend ) +
  geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
  scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
  facet_wrap( ~weight_class ) +
  geom_point(
    data = previsoes %>% filter(.fitted > 0.5), 
    aes(y = red_mais_velho, x = idade_red ), 
    color = "lightblue",
    alpha = 0.1
  ) +
  geom_point(
    data = previsoes %>% filter(.fitted < 0.5), 
    aes(y = red_mais_velho, x = idade_red ), 
    color = "lightcoral",
    alpha = 0.1
  ) +
  theme_minimal()

Rodando várias especificações de uma vez com `purrr`

Podemos treinar o modelo para cada classe usando group_by + nest() + map().

Estes comandos criam um tibble com uma coluna de tibbles aninhados já separados por classe de peso. A função map() pode, então rodar o modelo para cada classe.

modelo_ufc_por_classe <- ufc_aprend %>% 
  group_by(weight_class) %>% 
  nest_legacy() %>%
  mutate(
    modelo = map(data, ~lm(zebra_blue ~ idade_red  + red_mais_velho, data = .x)),
  ) %>% 
  select(-data)

ufc_por_classe <- modelo_ufc_por_classe %>%
  mutate(
    aumentado = map(modelo, augment)
  ) %>%
  unnest(aumentado)

ufc_por_classe_treinamento  <- ufc_por_classe %>%
  ungroup() %>%
  mutate(
    previsao_zebra_blue = as_factor(round(.fitted)),
    zebra_blue = as_factor(zebra_blue)
  )



confusionMatrix(ufc_por_classe_treinamento$previsao_zebra_blue, ufc_por_classe_treinamento$zebra_blue )

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 3128 1376
##          1  152  212
##                                           
##                Accuracy : 0.6861          
##                  95% CI : (0.6729, 0.6991)
##     No Information Rate : 0.6738          
##     P-Value [Acc > NIR] : 0.03414         
##                                           
##                   Kappa : 0.1088          
##                                           
##  Mcnemar's Test P-Value : < 2e-16         
##                                           
##             Sensitivity : 0.9537          
##             Specificity : 0.1335          
##          Pos Pred Value : 0.6945          
##          Neg Pred Value : 0.5824          
##              Prevalence : 0.6738          
##          Detection Rate : 0.6426          
##    Detection Prevalence : 0.9252          
##       Balanced Accuracy : 0.5436          
##                                           
##        'Positive' Class : 0               
##

Plotando as várias especificações de regressões lineares

grid_previsao_com_modelo <- grid_previsao %>%
  group_by(weight_class) %>% 
  nest() %>% 
  inner_join(modelo_ufc_por_classe, by = c("weight_class")) %>% 
  mutate(previsoes = map2(.x = modelo, .y =  data, .f =  ~augment(x = .x, newdata = .y) )) %>% 
  unnest_legacy(previsoes)



ggplot(ufc_aprend ) +
  geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
  scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
  facet_wrap( ~weight_class ) +
  geom_point(
    data = grid_previsao_com_modelo %>% filter(.fitted > 0.5), 
    aes(y = red_mais_velho, x = idade_red ), 
    color = "lightblue",
    alpha = 0.1
  ) +
  geom_point(
    data = grid_previsao_com_modelo %>% filter(.fitted < 0.5), 
    aes(y = red_mais_velho, x = idade_red ), 
    color = "lightcoral",
    alpha = 0.1
  ) +
  theme_minimal()

Separando dados de treinamento e validação

set.seed(2512)

ufc_aprend_classes_selec <- ufc_aprend %>% 
  filter(weight_class %in%
      c(   
      "Flyweight",
      "Bantamweight",
      "Featherweight",
      "Lightweight",
      "Welterweight",
      "Middleweight",
      "Light Heavyweight",
      "Heavyweight"
      )
      )



index_treino <- createDataPartition(
  ufc_aprend_classes_selec$zebra_blue,
  p = 0.75,
  list = FALSE,
  times = 1
) 


ufc_aprend_treino <- ufc_aprend_classes_selec %>% 
  filter(row_number() %in% index_treino) %>% 
  select(
    idade_red,
    red_mais_velho,
    zebra_blue,
    weight_class
  )

ufc_aprend_teste <- ufc_aprend_classes_selec %>% 
  filter(!(row_number() %in% index_treino)) %>% 
  select(
    idade_red,
    red_mais_velho,
    zebra_blue,
    weight_class
  )
  


modelo_ufc_por_classe_treino  <-  ufc_aprend_treino %>% 
  group_by(weight_class) %>% 
  nest_legacy() %>% 
  mutate( 
    modelo = map(data, ~lm(zebra_blue ~ idade_red  + red_mais_velho, data = .x)),
  ) %>% 
  select(-data)


resultados_teste_lm <-  ufc_aprend_teste %>% 
  group_by(weight_class) %>% 
  nest_legacy() %>% 
  inner_join(modelo_ufc_por_classe_treino, by = c("weight_class") ) %>% 
  mutate(
    previsao = map2(.x = modelo, .y =  data,  .f = ~predict.lm(.x, .y)  )
  ) %>%
  unnest(cols = c("data","previsao")) %>% 
  ungroup() %>% 
  mutate(
    actual = if_else(zebra_blue == 1, "B", "R"),
    predicted = if_else(previsao > 0.5, "B", "R"),
  ) %>% 
  mutate(
    actual  = fct_relevel(actual, "R", "B"),
    predicted = fct_relevel(predicted, "R", "B")
  ) %>% 
  select(
    actual,
    predicted
  ) 


confusionMatrix(resultados_teste_lm$predicted, resultados_teste_lm$actual  )

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   R   B
##          R 733 325
##          B  31  38
##                                           
##                Accuracy : 0.6841          
##                  95% CI : (0.6561, 0.7112)
##     No Information Rate : 0.6779          
##     P-Value [Acc > NIR] : 0.3405          
##                                           
##                   Kappa : 0.0814          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.9594          
##             Specificity : 0.1047          
##          Pos Pred Value : 0.6928          
##          Neg Pred Value : 0.5507          
##              Prevalence : 0.6779          
##          Detection Rate : 0.6504          
##    Detection Prevalence : 0.9388          
##       Balanced Accuracy : 0.5321          
##                                           
##        'Positive' Class : R               
##

Rodando um modelo mais complexo: Generalized Additive Model

#ALGO ERRADO AO RODAR O MARKDOWN. NÃO SEI



# modelo_ufc_por_classe_treino_gam  <-  ufc_aprend_treino %>%
#   # group_by(weight_class) %>%
#   nest(data =c(zebra_blue, idade_red, red_mais_velho ) ) %>%
#   mutate(
#     modelo = map(data, ~gam(formula = zebra_blue ~ s(idade_red) + s(red_mais_velho), data = .x))
#   ) %>%
#   select(-data)
# 
# 
# resultados_teste_gam <-  ufc_aprend_teste %>%
#   group_by(weight_class) %>%
#   nest_legacy() %>%
#   inner_join(modelo_ufc_por_classe_treino_gam, by = c("weight_class") ) %>%
#   mutate(
#     previsao = map2(.x = modelo, .y =  data,  .f = ~predict(.x, .y)  )
#   ) %>%
#   unnest(cols = c("data","previsao")) %>%
#   ungroup() %>%
#   mutate(
#     actual = if_else(zebra_blue == 1, "B", "R"),
#     predicted = if_else(previsao > 0.5, "B", "R"),
#   ) %>%
#   mutate(
#     actual  = fct_relevel(actual, "R", "B"),
#     predicted = fct_relevel(predicted, "R", "B")
#   ) %>%
#   select(
#     actual,
#     predicted
#   )
# 
# write_rds(resultados_teste_gam, "dados/ufc/resultados_teste_gam.rds")

resultados_teste_gam <- read_rds("dados/ufc/resultados_teste_gam.rds")

confusionMatrix(resultados_teste_gam$predicted, resultados_teste_gam$actual  )

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   R   B
##          R 716 310
##          B  48  53
##                                           
##                Accuracy : 0.6823          
##                  95% CI : (0.6543, 0.7095)
##     No Information Rate : 0.6779          
##     P-Value [Acc > NIR] : 0.3884          
##                                           
##                   Kappa : 0.1026          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.9372          
##             Specificity : 0.1460          
##          Pos Pred Value : 0.6979          
##          Neg Pred Value : 0.5248          
##              Prevalence : 0.6779          
##          Detection Rate : 0.6353          
##    Detection Prevalence : 0.9104          
##       Balanced Accuracy : 0.5416          
##                                           
##        'Positive' Class : R               
##

Plotando as regiões deste modelo mais complexo

É possível observar que este modelo se adapta com muito mais potência às especificidades dos dados. Ou seja, consegue um viés baixo no conjunto de treinamento.

Ele se adapta tão bem ao conjunto que se apresentarmos um conjunto um pouco diferente, a \(\hat{f}()\) será diferente certamente.

O modelo parece se adaptar demais aos dados, gerando overfitting. Será que ele é melhor do que um modelo simples?

Precisamos de um mecanismo para avaliar o poder de generalização do modelo.

# 
# grid_previsao_com_modelo_gam <- grid_previsao %>%
#   group_by(weight_class) %>% 
#   nest() %>% 
#   inner_join(modelo_ufc_por_classe_treino_gam, by = c("weight_class")) %>% 
#   mutate(previsoes = map2(.x = modelo, .y =  data, .f =  ~predict( .x, .y) )) %>% 
#   select(-modelo) %>% 
#   unnest(cols = c("previsoes","data")) %>% 
#   ungroup()
# 
# 
# write_rds(grid_previsao_com_modelo_gam, "dados/ufc/grid_previsao_com_modelo_gam.RDS") 

grid_previsao_com_modelo_gam <- read_rds("dados/ufc/grid_previsao_com_modelo_gam.RDS")



ggplot( bind_rows(ufc_aprend_classes_selec) ) +
  geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
  scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
  facet_wrap( ~weight_class ) +
  geom_point(
    data = grid_previsao_com_modelo_gam %>% filter(previsoes > 0.5),
    aes(y = red_mais_velho, x = idade_red ),
    color = "lightblue",
    alpha = 0.1
  ) +
  geom_point(
    data = grid_previsao_com_modelo_gam %>% filter(previsoes < 0.5),
    aes(y = red_mais_velho, x = idade_red ),
    color = "lightcoral",
    alpha = 0.1
  ) +
  theme_minimal()

Execução de modelos com caret

Vamos mostrar a biblioteca caret, que facilita o processo de aprendizado de vários modelos com várias especificações de hiperparâmetros, executando o processo de cross-validation.

Primeiro vamos ler uma base com informações e diagnósticos de câncer de mama (UCI 2019)

Já vamos separar os dados de teste

options(OutDec = ",")


dados <- read_csv("dados/diagnostico/wdbc.csv") %>% 
    select(-"ID number") %>% 
    rename_all( .funs = ~str_replace(.,"fractal-media ","fractal_")   ) %>% 
    rename_all( .funs = ~str_replace(.,"fractal-pior ","fractal_")   ) %>% 
    rename_all( .funs = ~str_replace(.,"fractal-dv ","fractal_")   ) %>% 
    rename("fractal_dimension-media" = "fractal_dimension-meia") %>% 
    rename_all( .funs = ~str_replace(.,"-","_")   ) %>% 
    rename_all( .funs = ~str_replace(.," ","_")   ) %>% 
    mutate(Diagnosis = as.factor(Diagnosis))

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Diagnosis = col_character()
## )

## See spec(...) for full column specifications.

set.seed(13)



rows <- sample(nrow(dados))

dados <- dados[rows,]

split <- round(nrow(dados) * .70 )

treino <- dados[1:split, ]

teste <- dados[(split + 1):nrow(dados), ]

Preparando 4-fold cross validation

A função trainControl() estabelece os conjuntos para cross validation.

Os mesmos conjuntos serão usados para todos os modelos

set.seed(13)


controle_cv <- trainControl(
    
    summaryFunction = twoClassSummary,
    classProbs = TRUE,
    verboseIter = FALSE,
    returnResamp = "all",
    number = 4,
    repeats = 10,
    returnData = TRUE,
    savePredictions = "all",
    method = "repeatedcv",
    allowParallel = FALSE,

)

Treinando um modelo de classificação. Exemplo: Regressão logística

A função train() permite o estabelecimento de uma métrica para a escolha do melhor valor para os hiperparâmetros (no caso, não existem), e execuções de funções de pré-processamento.

A função retorna um objeto com todas as informações relativas ao treinamento e a melhor especificação do modelo treinada com todos os dados passados.

A saída padrão impressa em problemas de classificação mostra o resultado do melhor modelo nos seus dados de validação dele.

model_logistic <- train( 
    
    form = Diagnosis ~ . ,
    data = treino,
    metric = "ROC",
    method = "glm",
    trControl = controle_cv,
    preProcess = c( "center", "scale"),
    family=binomial(link='logit')

)


model_logistic

## Generalized Linear Model 
## 
## 398 samples
##  30 predictor
##   2 classes: 'B', 'M' 
## 
## Pre-processing: centered (30), scaled (30) 
## Resampling: Cross-Validated (4 fold, repeated 10 times) 
## Summary of sample sizes: 299, 298, 299, 298, 299, 298, ... 
## Resampling results:
## 
##   ROC        Sens       Spec     
##   0,9549433  0,9515084  0,9121429

Avaliando a curva ROC

A curva ROC (Receiver Operating Characteristics) foi inventada na época da Segunda Guerra Mundial para avaliar se os operadores de radar americanos estavam detectando confiavelmente aeronaves japonesas a partir de sinais de radar.

A curva mostra, para vários thresholds, qual a fração de verdadeiros positivos (ou sensibilidade) e a fração de falsos positivos (fall-out, ou \(1 - especificidade\) ).

Uma métrica numérica que traduz o quão bom um modelo de classificação é consiste na área embaixo desta curva (AUC, Area Under the Curve). Note que quanto mais perto de um essa área, menor a taxa de falsos positivos e maior a sensibilidade

ggplot(model_logistic$pred, aes(m = M, d = obs, color = Resample )) +
    geom_roc( labels = FALSE  ) +
    coord_equal() + style_roc() + ggtitle("ROC", subtitle = "Métricas para diversos thresholds" )+       style_roc()

Desempenho dos diversos modelos

Podemos avaliar a ROC de cada treinamento feito.

result_model_logistic <-  model_logistic$resample %>% 
    select(Resample, ROC) %>%
    rename(AUC = ROC) 


result_model_logistic %>% 
    mutate_if(is.numeric, percent) %>% 
    kable( caption = "\\label{tab_auc_reglog}Métricas para cada Fold da regressão logistica")

Métricas para cada Fold da regressão logistica
Resample	AUC
Fold1.Rep01	94.4%
Fold2.Rep01	96.1%
Fold3.Rep01	93.0%
Fold4.Rep01	95.6%
Fold1.Rep02	96.5%
Fold2.Rep02	93.8%
Fold3.Rep02	95.1%
Fold4.Rep02	95.3%
Fold1.Rep03	95.3%
Fold2.Rep03	96.7%
Fold3.Rep03	94.2%
Fold4.Rep03	96.7%
Fold1.Rep04	94.2%
Fold2.Rep04	95.4%
Fold3.Rep04	95.9%
Fold4.Rep04	96.5%
Fold1.Rep05	98.1%
Fold2.Rep05	99.2%
Fold3.Rep05	92.0%
Fold4.Rep05	96.4%
Fold1.Rep06	98.4%
Fold2.Rep06	90.0%
Fold3.Rep06	97.1%
Fold4.Rep06	96.9%
Fold1.Rep07	95.1%
Fold2.Rep07	98.1%
Fold3.Rep07	96.8%
Fold4.Rep07	88.0%
Fold1.Rep08	96.0%
Fold2.Rep08	90.7%
Fold3.Rep08	96.6%
Fold4.Rep08	95.4%
Fold1.Rep09	96.5%
Fold2.Rep09	95.2%
Fold3.Rep09	95.3%
Fold4.Rep09	96.1%
Fold1.Rep10	95.2%
Fold2.Rep10	98.4%
Fold3.Rep10	95.7%
Fold4.Rep10	97.7%

Desempenho e estimativa da variância do modelo

result_model_logistic %>%     
    summarise(media = mean(AUC), sd = sd(AUC)) %>% 
    rename("Média AUC" = media, "Desvio-padrão AUC" = sd) %>% 
    mutate_if(is.numeric, percent) %>% 
    kable(caption = "\\label{tab_auc_reglog_grup}Média e desvio-padrão da Métrica AUC para regressão logistica")

Média e desvio-padrão da Métrica AUC para regressão logistica
Média AUC	Desvio-padrão AUC
95.5%	2.27%

Gerando um latex bonitinho do modelo

texreg(model_logistic$finalModel, custom.model.names = c("Regressão logística simples"), caption = "Coeficientes estimados na regressão logística", label = "reg:log", fontsize = "footnotesize")

## 
## \begin{table}
## \begin{center}
## \begin{footnotesize}
## \begin{tabular}{l c }
## \hline
##  & Regressão logística simples \\
## \hline
## (Intercept)               & $80,96$        \\
##                           & $(35914,01)$   \\
## radius\_media             & $-264,43$      \\
##                           & $(767230,04)$  \\
## texture\_media            & $41,40$        \\
##                           & $(21122,91)$   \\
## perimeter\_media          & $1183,18$      \\
##                           & $(1099949,43)$ \\
## area\_media               & $-1040,37$     \\
##                           & $(494524,57)$  \\
## smoothness\_media         & $-65,13$       \\
##                           & $(30919,21)$   \\
## compactness\_media        & $-364,13$      \\
##                           & $(79670,45)$   \\
## concavity\_media          & $6,25$         \\
##                           & $(189198,65)$  \\
## concave\_points\_media    & $367,79$       \\
##                           & $(89897,83)$   \\
## symmetry\_media           & $26,33$        \\
##                           & $(26581,88)$   \\
## fractal\_dimension\_media & $80,71$        \\
##                           & $(36331,35)$   \\
## radius\_dv                & $246,29$       \\
##                           & $(104667,71)$  \\
## texture\_dv               & $14,71$        \\
##                           & $(31824,49)$   \\
## perimeter\_dv             & $73,35$        \\
##                           & $(67727,22)$   \\
## area\_dv                  & $-468,01$      \\
##                           & $(232868,82)$  \\
## smoothness\_dv            & $-73,73$       \\
##                           & $(23297,51)$   \\
## compactness\_dv           & $100,76$       \\
##                           & $(46184,65)$   \\
## concavity\_dv             & $-24,23$       \\
##                           & $(67029,28)$   \\
## concave\_points\_dv       & $51,13$        \\
##                           & $(46579,32)$   \\
## symmetry\_dv              & $71,34$        \\
##                           & $(24137,24)$   \\
## fractal\_dimension\_dv    & $-236,13$      \\
##                           & $(68224,41)$   \\
## radius\_pior              & $-817,84$      \\
##                           & $(427324,69)$  \\
## texture\_pior             & $0,44$         \\
##                           & $(18864,24)$   \\
## perimeter\_pior           & $-187,19$      \\
##                           & $(129752,63)$  \\
## area\_pior                & $1626,90$      \\
##                           & $(633836,01)$  \\
## smoothness\_pior          & $56,67$        \\
##                           & $(25600,89)$   \\
## compactness\_pior         & $-114,58$      \\
##                           & $(56734,42)$   \\
## concavity\_pior           & $95,60$        \\
##                           & $(74790,04)$   \\
## concave\_points\_pior     & $-22,85$       \\
##                           & $(41737,35)$   \\
## symmetry\_pior            & $-28,93$       \\
##                           & $(38034,29)$   \\
## fractal\_dimension\_pior  & $203,37$       \\
##                           & $(76158,01)$   \\
## \hline
## AIC                       & 62,00          \\
## BIC                       & 185,58         \\
## Log Likelihood            & -0,00          \\
## Deviance                  & 0,00           \\
## Num. obs.                 & 398            \\
## \hline
## \multicolumn{2}{l}{\tiny{$^{***}p<0,001$, $^{**}p<0,01$, $^*p<0,05$}}
## \end{tabular}
## \end{footnotesize}
## \caption{Coeficientes estimados na regressão logística}
## \label{reg:log}
## \end{center}
## \end{table}

Treinando outro modelo: lasso-ridge

Uma outra forma de evitar que muitas variáveis sejam usadas no modelo é aplicar uma penalidade diminuir a variância do modelo. É isso que as regressões do tipo Ridge e Lasso fazem.

O modelo mostrado nessa seção é rodado com várias parametrizações. Este modelo, chamado Elastic Net, conjuga a penalização do tipo Ridge com a penalização do tipo Lasso modificando a função de penalização da regressão, que originalmente é o erro quadrático:

\[RSS = \sum_{i = 1}^{n} ( y_i - \beta_0 - \sum_{j=1}^{p}\beta_j x_{ij})^2 \]

Para a regressão Ridge, os coeficientes são penalizados de forma quadrática. Isso diminui a variância do modelo mas não diminui tanto a diminuição do número de coeficientes diferentes de 0:

\[Loss_{Ridge} = RSS + \lambda \sum_{j=1}^{p}\beta_j^2 \]

Para a regressão Lasso, os coeficientes são penalizados pelo seu valor absoluto. Isso diminui a variância do modelo E diminui o número de coeficientes diferentes de 0, favorecendo a interpretabilidade

\[Loss_{Lasso} = RSS + \lambda \sum_{j=1}^{p} \left| \beta_j \right| \]

Conjugando Lasso e Ridge e avaliando o modelo

Cada conjunto de hiperparâmetros do modelo Elastic Net rodado nesta possui dois valores. Um dos valores, \(\alpha\) define que peso deve ser dado a cada tipo de penalização Ridge ou Lasso, onde \(\alpha = 0\) é Ridge puro e \(\alpha = 1\) é Lasso puro. O outro parâmetro, \(\lambda\), regula a intensidade da penalização.

Os resultados para várias parametrizações é mostrado abaixo.

Veja como podemos passar um grid de hiperparâmetros e avaliar esses resultados.

set.seed(13)

model_net <- train(
    
    Diagnosis ~ .,
    treino,
    metric = "ROC",
    method = "glmnet",
    trControl = controle_cv,
    preProcess = c( "center", "scale"),
    tuneGrid = expand.grid(
         alpha = c(0,0.5,0.75,0.1,0.125,0.15,0.25,0.5, 1),
         lambda = 0:10/500
     )    
    

)


plot(model_net)

Os resultados de cada parametrização são mostrados em ordem decrescente de AUC.

resultados_net <- model_net$resample %>% 
    group_by(alpha, lambda) %>% 
    summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>% 
    arrange(desc(media)) %>% 
    rename("Média AUC" = media) 


head(resultados_net, 20)

## # A tibble: 20 x 4
## # Groups:   alpha [6]
##    alpha lambda `Média AUC` `Desvio-padrão AUC`
##    <dbl>  <dbl>       <dbl>               <dbl>
##  1 1      0.02        0.994             0.00746
##  2 1      0.014       0.994             0.00762
##  3 1      0.018       0.994             0.00756
##  4 1      0.016       0.994             0.00761
##  5 1      0.012       0.994             0.00776
##  6 1      0.01        0.994             0.00782
##  7 1      0.008       0.994             0.00777
##  8 0.75   0.006       0.994             0.00790
##  9 0.1    0.01        0.994             0.00715
## 10 0.25   0.01        0.994             0.00732
## 11 0.125  0.018       0.994             0.00732
## 12 0.15   0.008       0.994             0.00720
## 13 0.1    0.02        0.994             0.00730
## 14 1      0.006       0.994             0.00784
## 15 0.1    0.016       0.994             0.00724
## 16 0.1    0.018       0.994             0.00726
## 17 0.15   0.014       0.994             0.00728
## 18 0.1    0.012       0.994             0.00716
## 19 0.75   0.004       0.994             0.00766
## 20 0.125  0.008       0.994             0.00721

Pré-processando com Principal Component Analysis

Uma outra forma de reduzir a dimensionalidade do probelam é usar o método PCA (Principal Component Analysis) para transformar as variáveis em componentes principais em menor número mas com boa parte da informação original.

Veja como estabelecemos esse pré-processamento.

set.seed(13)

model_net_pca <- train(
    
    Diagnosis ~ .,
    treino,
    metric = "ROC",
    method = "glmnet",
    trControl = controle_cv,
    preProcess = c( "center", "scale", "pca"),
    tuneGrid = expand.grid(
         alpha = c(0,0.5,0.75,0.1,0.125,0.15,0.25,0.5, 1),
         lambda = 0:10/500
     )    
    

)


plot(model_net_pca)

resultados_net_pca <- model_net_pca$resample %>% 
    group_by(alpha, lambda) %>% 
    summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>% 
    arrange(desc(media)) %>% 
    rename("Média AUC" = media)

resultados_net_pca

## # A tibble: 88 x 4
## # Groups:   alpha [8]
##    alpha lambda `Média AUC` `Desvio-padrão AUC`
##    <dbl>  <dbl>       <dbl>               <dbl>
##  1  1     0.01        0.993             0.00664
##  2  1     0.012       0.993             0.00658
##  3  1     0.014       0.993             0.00648
##  4  1     0.008       0.993             0.00692
##  5  0.75  0.016       0.993             0.00709
##  6  0.75  0.012       0.993             0.00734
##  7  0.75  0.014       0.993             0.00719
##  8  0.75  0.006       0.993             0.00754
##  9  0.75  0.008       0.993             0.00758
## 10  1     0.006       0.993             0.00703
## # ... with 78 more rows

Árvore de decisão

A árvore de decisão particiona o espaço formado pelas variáveis explicativas em subespaços baseando-se na “pureza” desses subespaços com relação à variável dependente.

model_tree <- train(
    
    Diagnosis  ~ .,
    treino,
    metric = "ROC",
    method = "rpart",
    trControl = controle_cv,
    preProcess = c( "center", "scale")
    


)


resultados_tree <- model_tree$resample %>% 
    group_by(cp) %>% 
    summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>% 
    arrange(desc(media)) %>% 
    rename("Média AUC" = media)


resultados_tree

## # A tibble: 3 x 3
##       cp `Média AUC` `Desvio-padrão AUC`
##    <dbl>       <dbl>               <dbl>
## 1 0.0214       0.942              0.0273
## 2 0.0500       0.932              0.0293
## 3 0.829        0.773              0.193

Random Forests

A forma como a árvore de decisão é criada faz com que ela tenha muita variância.

Cada decisão de particionamento é tomada a partir de características que podem ser muito específicas aos dados de treinamento.

Uma ideia usada nas Random Forests é criar conjuntos de treinamento criados a partir do conjunto original, mas retirando amostras de mesmo tamanho COM REPOSIÇÃO. Além disso, cada vez que a árvore é particionada, a partição só pdoe acontecer em \(m\) das variáveis explicativas.

O resultado final é uma média do resultado dessas várias árvores

Estas duas mudanças fazem com que o modelo tenha uma variância muito menor do que as árvore de decisão simples

model_ranger <- train(
    
    Diagnosis  ~ .,
    treino,
    metric = "ROC",
    method = "ranger",
    trControl = controle_cv,
    preProcess = c( "center", "scale"),
    verbose = TRUE,
    tuneGrid = expand.grid(
         mtry = seq(2, by = 2, to = 20),
         splitrule = c("gini", "extratrees"),
         min.node.size = 1
     )    


)

plot(model_ranger)

resultados_ranger <- model_ranger$resample %>% 
    group_by(mtry, splitrule) %>% 
    summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>% 
    arrange(desc(media)) %>% 
    rename("Média AUC" = media)


resultados_ranger

## # A tibble: 20 x 4
## # Groups:   mtry [10]
##     mtry splitrule  `Média AUC` `Desvio-padrão AUC`
##    <dbl> <fct>            <dbl>               <dbl>
##  1     2 extratrees       0.992             0.00606
##  2     4 extratrees       0.992             0.00647
##  3     8 extratrees       0.992             0.00682
##  4     6 extratrees       0.991             0.00686
##  5    16 extratrees       0.991             0.00778
##  6    10 extratrees       0.991             0.00759
##  7    12 extratrees       0.990             0.00808
##  8    14 extratrees       0.990             0.00847
##  9    20 extratrees       0.990             0.00848
## 10    18 extratrees       0.990             0.00913
## 11     2 gini             0.989             0.00933
## 12     6 gini             0.988             0.0104 
## 13    18 gini             0.988             0.0101 
## 14    16 gini             0.988             0.0101 
## 15    10 gini             0.988             0.0104 
## 16     4 gini             0.988             0.0105 
## 17    20 gini             0.988             0.0101 
## 18     8 gini             0.987             0.0106 
## 19    14 gini             0.987             0.0104 
## 20    12 gini             0.987             0.0108

Montando um modelo ensemble

Podemos montar um modelo ensemble, onde cada modelo original vota e a escolha final recai sobre a classificação que receber mais votos

pred_logistic <- model_logistic$pred %>% 
    semi_join(model_logistic$bestTune) %>% 
    mutate(modelo = "Logística")

pred_net_pca <- model_net_pca$pred %>% 
    semi_join(model_net_pca$bestTune) %>% 
    mutate(modelo = "Lasso-Ridge PCA")
    
pred_ranger <- model_ranger$pred %>% 
    semi_join(model_ranger$bestTune) %>% 
    mutate(modelo = "Floresta")

pred_ensemble_todos <- bind_rows(pred_logistic, pred_net_pca, pred_ranger)

pred_ensemble_mediana <- pred_ensemble_todos %>% 
    separate(Resample, c("Fold", "Repeat")) %>% 
    group_by(Repeat, rowIndex) %>% 
    summarise(M = median(M), B= median(B), n = n(), obs = unique(obs)) 

ret_roc <- function(data)
{
    roc(data$obs,data$M)
}


ret_roc_auc <- function(data)
{
    unlist(as.numeric(roc(data$obs,data$M)$auc))
}


pred_ensemble_mediana_foldado <- model_ranger$pred %>% 
    semi_join(model_ranger$bestTune) %>% 
    separate(Resample, c("Fold", "Repeat")) %>% 
    select(Fold, Repeat, rowIndex) %>% 
    inner_join(pred_ensemble_mediana, by = c("Repeat", "rowIndex"), suffix = c("","_") )  %>% 
    group_by(Fold, Repeat) %>% 
    nest() %>% 
    mutate( roc = map(data, ret_roc), auc = unlist(map(data, ret_roc_auc))  )

resultado_ensemble_pra_mostrar <- tibble( "Média AUC" = mean(pred_ensemble_mediana_foldado$auc), "Desvio-padrão AUC" = sd(pred_ensemble_mediana_foldado$auc) )


resultado_ensemble_pra_mostrar

## # A tibble: 1 x 2
##   `Média AUC` `Desvio-padrão AUC`
##         <dbl>               <dbl>
## 1       0.992             0.00748

Uma forma de avaliar o impacto das variáveis

É possível propor uma forma de avaliação do impacto das variáveis indireta e agnóstica ao modelo, ou seja, que pode ser usada da mesma forma para qualquer modelo e .

Para avaliar o impacto de cada variável nos resultados de cada um dos 5 modelos, podemos gerar casos sintéticos e avaliar a probabilidade de esses casos serem malignos segundo cada modelo.

Dois grupos de casos sintéticos foram gerados.

No primeiro grupo, para cada variável explicativa \(x_i\), foram criados casos onde o valor de \(x_i\) varia em relação a seus quantis, enquanto as outras permanecem parada no valor das suas medianas.

No segundo grupo, para avaliar uma região onde os casos têm mais chance de ser malignos, as variáveis que não estão sendo avaliadas a cada gráfico são mantidas no quantil 0,65.

medias_treino <- treino %>% 
    select(-Diagnosis) %>% 
    summarise_all(mean) %>% 
    gather(variavel, media) 

dps_treino <- treino %>% 
    select(-Diagnosis) %>% 
    summarise_all(sd) %>% 
    gather(variavel, dp) 

mediana_treino <- treino %>% 
    select(-Diagnosis) %>% 
    summarise_all(median) %>% 
    gather(variavel, mediana) 


variaveis <- medias_treino %>% 
    select(variavel)

caso_media <- medias_treino %>% 
    spread(variavel, media)


quantis <- treino %>% 
    select(-Diagnosis) %>% 
    gather(variavel, valor) %>% 
    crossing( tibble( q = seq(0,1, 0.05))   ) %>% 
    group_by(variavel,q ) %>% 
    summarise( valor = quantile(valor, mean(q))) 
    

casos <- variaveis %>% 
    crossing(variaveis %>% rename(fixa = variavel)) %>% 
    crossing(tibble( q = seq(0,1, 0.05)) ) %>% 
    arrange(q, variavel, fixa) %>% 
    mutate(caso = (row_number()-1) %/% nrow(variaveis) + 1) %>% 
    mutate(q_cada_variavel = if_else(variavel == fixa, q, 0.5 )) %>% 
    inner_join(quantis, by=c("q_cada_variavel" = "q", "fixa" = "variavel") ) %>%
    select(-q_cada_variavel) %>% 
    spread(fixa, valor)
    
prev_log <- predict(model_logistic, casos, type = "prob") %>% 
    mutate(modelo = "Logistica") %>% 
    bind_cols(casos)

prev_ranger <- predict(model_ranger, casos, type = "prob") %>% 
    mutate(modelo = "Floresta") %>% 
    bind_cols(casos)

prev_pca <- predict(model_net_pca, casos, type = "prob") %>% 
    mutate(modelo = "Lasso-Ridge PCA") %>% 
    bind_cols(casos)



prev <- 
    bind_rows(prev_log, prev_ranger, prev_pca)



prev %>% 
    ggplot() +
    geom_line( aes(x = q, y = M, color = modelo) ) +
    facet_wrap( ~variavel) +
    theme_minimal() +
    theme(
        aspect.ratio = .4,
        strip.text = element_text(size = 5),
        strip.background.y = element_rect(size = c(0.1,0.1)),
        strip.switch.pad.grid = unit(.01,"cm"),
        strip.switch.pad.wrap = unit(.01,"cm"),
        axis.text =  element_text(size = 4) 
        
        ) +
  theme(legend.position = "top" )

casos <- variaveis %>% 
    crossing(variaveis %>% rename(fixa = variavel)) %>% 
    crossing(tibble( q = seq(0,1, 0.05)) ) %>% 
    arrange(q, variavel, fixa) %>% 
    mutate(caso = (row_number()-1) %/% nrow(variaveis) + 1) %>% 
    mutate(q_cada_variavel = if_else(variavel == fixa, q, 0.65 )) %>% 
    inner_join(quantis, by=c("q_cada_variavel" = "q", "fixa" = "variavel") ) %>%
    select(-q_cada_variavel) %>% 
    spread(fixa, valor)
    
prev_log <- predict(model_logistic, casos, type = "prob") %>% 
    mutate(modelo = "Logistica") %>% 
    bind_cols(casos)

prev_ranger <- predict(model_ranger, casos, type = "prob") %>% 
    mutate(modelo = "Floresta") %>% 
    bind_cols(casos)

prev_pca <- predict(model_net_pca, casos, type = "prob") %>% 
    mutate(modelo = "Lasso-Ridge PCA") %>% 
    bind_cols(casos)



prev <- 
    bind_rows(prev_log, prev_ranger, prev_pca)



prev %>% 
    ggplot() +
    geom_line( aes(x = q, y = M, color = modelo) ) +
    facet_wrap( ~variavel) +
    theme_minimal() +
    theme(
        aspect.ratio = .4,
        strip.text = element_text(size = 5),
        strip.background.y = element_rect(size = c(0.1,0.1)),
        strip.switch.pad.grid = unit(.01,"cm"),
        strip.switch.pad.wrap = unit(.01,"cm"),
        axis.text =  element_text(size = 4) 
        
        ) +
  theme(legend.position = "top" )

Modelos escolhidos e seus thresholds: Model Selection

roc_best_modelo <- function(modelo)
{

    dados_roc <- modelo$pred %>% 
        semi_join(modelo$bestTune) %>% 
        mutate(obs = if_else(obs == "B", 0, 1))

    dados_curva <- roc(dados_roc$obs, dados_roc$M)
    
    tibble( sensitivities = dados_curva$sensitivities, 
            specificities= dados_curva$specificities,  
            thresholds = dados_curva$thresholds,
            auc = dados_curva$auc
            ) 
    
}   


pred_best_modelo <- function(modelo)
{
    modelo$pred %>% 
        semi_join(modelo$bestTune)

}

    
plota_roc_modelo <-  function(modelo, prob, nome)
{
    
    
    
    
    roc_modelo <-  roc_best_modelo(modelo)
    pred_modelo <- pred_best_modelo(modelo)
    
    
    threshold_escolhido <- roc_modelo %>% 
        filter(thresholds >= prob) %>% 
        top_n(n = 1, wt = desc(thresholds )) 
    
    
    ggplot(pred_modelo ) +
        geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs )  ) +
        geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
        geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
        coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
    
}


plota_roc_ensemble <-  function(modelo, prob, nome)
{
    prob <- 0.3
    nome <- "Ensemble"
    dados_curva <- roc(pred_ensemble_mediana$obs, pred_ensemble_mediana$M)
    
    roc_modelo <- tibble( sensitivities = dados_curva$sensitivities, 
            specificities= dados_curva$specificities,  
            thresholds = dados_curva$thresholds,
            auc = dados_curva$auc
            ) 
    
    
    
    pred_modelo <- pred_ensemble_mediana
    
    
    threshold_escolhido <- roc_modelo %>% 
        filter(thresholds >= prob) %>% 
        top_n(n = 1, wt = desc(thresholds )) 
    
    
    ggplot(pred_modelo ) +
        geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs )  ) +
        geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
        geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
        coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
    
}


plota_roc_modelo(model_logistic, 1e-8, "Logistica")

plota_roc_modelo(model_net_pca, 0.18, "Lasso-Ridge PCA")

plota_roc_modelo(model_ranger, 0.25, "Random Forest")

plota_roc_ensemble <-  function(modelo, prob, nome)
{
    prob <- 0.17
    nome <- "Ensemble"
    dados_curva <- roc(pred_ensemble_mediana$obs, pred_ensemble_mediana$M)
    
    roc_modelo <- tibble( sensitivities = dados_curva$sensitivities, 
            specificities= dados_curva$specificities,  
            thresholds = dados_curva$thresholds,
            auc = dados_curva$auc
            ) 
    
    
    
    pred_modelo <- pred_ensemble_mediana
    
    
    threshold_escolhido <- roc_modelo %>% 
        filter(thresholds >= prob) %>% 
        top_n(n = 1, wt = desc(thresholds )) 
    
    
    ggplot(pred_modelo ) +
        geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs )  ) +
        geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
        geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
        coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
    
}


plota_roc_ensemble()

Avaliação nos dados de teste: Model Assessment

confusao_teste <- function(modelo, threshold)
{

    prev_teste_svm <- predict(modelo, teste, type = "prob") 
    
    
    pred_teste_svm <- prev_teste_svm %>% 
        mutate(pred = as.factor( if_else(M > threshold, "M", "B"))) %>% 
        select(pred)
    
    conf_svm <- confusionMatrix(data = pred_teste_svm$pred, reference =  teste$Diagnosis, positive = "M")
}


gera_dados_confusao_teste <- function(modelo, threshold)
{

    prev_teste_svm <- predict(modelo, teste, type = "prob") 
    
    
    pred_teste_svm <- prev_teste_svm %>% 
        mutate(pred = as.factor( if_else(M > threshold, "M", "B"))) %>% 
        select(pred)
    
    #tibble(data = pred_teste_svm$pred, reference =  teste$Diagnosis, positive = "M")
    
    bind_cols(pred_teste_svm, teste) 
    
}

confusao_logistic <-  confusao_teste(model_logistic, 1e-8)

confusao_logistic

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 89  7
##          M 10 65
##                                          
##                Accuracy : 0,9006         
##                  95% CI : (0,8456, 0,941)
##     No Information Rate : 0,5789         
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0,7972         
##                                          
##  Mcnemar's Test P-Value : 0,6276         
##                                          
##             Sensitivity : 0,9028         
##             Specificity : 0,8990         
##          Pos Pred Value : 0,8667         
##          Neg Pred Value : 0,9271         
##              Prevalence : 0,4211         
##          Detection Rate : 0,3801         
##    Detection Prevalence : 0,4386         
##       Balanced Accuracy : 0,9009         
##                                          
##        'Positive' Class : M              
##

confusao_net_pca <-  confusao_teste(model_net_pca, 0.18)

confusao_net_pca

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 91  1
##          M  8 71
##                                           
##                Accuracy : 0,9474          
##                  95% CI : (0,9024, 0,9757)
##     No Information Rate : 0,5789          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0,8935          
##                                           
##  Mcnemar's Test P-Value : 0,0455          
##                                           
##             Sensitivity : 0,9861          
##             Specificity : 0,9192          
##          Pos Pred Value : 0,8987          
##          Neg Pred Value : 0,9891          
##              Prevalence : 0,4211          
##          Detection Rate : 0,4152          
##    Detection Prevalence : 0,4620          
##       Balanced Accuracy : 0,9527          
##                                           
##        'Positive' Class : M               
##

confusao_ranger  <-  confusao_teste(model_ranger, 0.25)

confusao_ranger

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 91  2
##          M  8 70
##                                           
##                Accuracy : 0,9415          
##                  95% CI : (0,8951, 0,9716)
##     No Information Rate : 0,5789          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0,8814          
##                                           
##  Mcnemar's Test P-Value : 0,1138          
##                                           
##             Sensitivity : 0,9722          
##             Specificity : 0,9192          
##          Pos Pred Value : 0,8974          
##          Neg Pred Value : 0,9785          
##              Prevalence : 0,4211          
##          Detection Rate : 0,4094          
##    Detection Prevalence : 0,4561          
##       Balanced Accuracy : 0,9457          
##                                           
##        'Positive' Class : M               
##

pred_test_1 <- predict(model_logistic, teste, type = "prob") %>% 
    mutate(row = row_number())  

pred_test_2 <- predict(model_net_pca, teste, type = "prob") %>% 
    mutate(row = row_number())  

pred_test_3 <- predict(model_ranger, teste, type = "prob") %>% 
    mutate(row = row_number())  

pred_teste_ensemble <- bind_rows(pred_test_1, pred_test_2, pred_test_3) %>% 
    group_by(row) %>% 
    summarise (M = median(M)) %>% 
    mutate(pred = as.factor( if_else(M > 0.17, "M", "B"))) 



confusao_ensemble <- confusionMatrix(data = pred_teste_ensemble$pred, reference =  teste$Diagnosis, positive = "M")
    
confusao_ensemble

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 91  1
##          M  8 71
##                                           
##                Accuracy : 0,9474          
##                  95% CI : (0,9024, 0,9757)
##     No Information Rate : 0,5789          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0,8935          
##                                           
##  Mcnemar's Test P-Value : 0,0455          
##                                           
##             Sensitivity : 0,9861          
##             Specificity : 0,9192          
##          Pos Pred Value : 0,8987          
##          Neg Pred Value : 0,9891          
##              Prevalence : 0,4211          
##          Detection Rate : 0,4152          
##    Detection Prevalence : 0,4620          
##       Balanced Accuracy : 0,9527          
##                                           
##        'Positive' Class : M               
##

dados_confusao_teste1 <- gera_dados_confusao_teste(model_logistic, 1e-8) %>% 
    mutate(modelo = "Logistic", linha = row_number())

dados_confusao_teste2 <- gera_dados_confusao_teste(model_net_pca, 0.18) %>% 
    mutate(modelo = "Lasso-Ridge", linha = row_number())

dados_confusao_teste3 <- gera_dados_confusao_teste(model_ranger , 0.25) %>% 
    mutate(modelo = "Random Forest", linha = row_number())


dados_confusao_teste_todos <- bind_rows(
    
    dados_confusao_teste1,
    dados_confusao_teste2,
    dados_confusao_teste3

)

pred_errados_teste <- dados_confusao_teste_todos %>% 
    filter(pred != Diagnosis) %>% 
    group_by(linha, pred, Diagnosis) %>% 
    summarise(n= n()) %>% 
    arrange(desc(n))

pred_errados_teste

## # A tibble: 23 x 4
## # Groups:   linha, pred [23]
##    linha pred  Diagnosis     n
##    <int> <fct> <fct>     <int>
##  1    37 M     B             3
##  2    47 M     B             3
##  3   111 B     M             3
##  4   148 M     B             3
##  5   150 M     B             3
##  6    94 M     B             2
##  7   120 M     B             2
##  8   139 M     B             2
##  9     2 B     M             1
## 10    10 M     B             1
## # ... with 13 more rows

Avaliando regiões de fraqueza do modelo

dados_confusao_teste_todos_tidy <- dados_confusao_teste_todos %>% 
    gather(atributo, valor, - pred,-Diagnosis, - modelo, - linha)

dados_confusao_teste_todos_certos <- 
    dados_confusao_teste_todos_tidy %>% 
    filter(pred == Diagnosis)

dados_confusao_teste_todos_falso_negativo <- 
    dados_confusao_teste_todos_tidy %>% 
    filter(pred == "M" & Diagnosis == "B")

dados_confusao_teste_todos_falso_positivo <- 
    dados_confusao_teste_todos_tidy %>% 
    filter(pred == "B" & Diagnosis == "M")


ggplot() +
    geom_jitter(data = dados_confusao_teste_todos_certos %>% filter(Diagnosis == "M"), aes(x = valor, y = 0), alpha = 0.05, size = 0.5, color = "red") +
    geom_jitter(data = dados_confusao_teste_todos_certos %>% filter(Diagnosis == "B")  , aes(x = valor, y = 0), alpha = 0.05, size = 0.5, color = "blue") +
    geom_jitter(data = dados_confusao_teste_todos_falso_negativo, aes(x = valor, y = 0), alpha = 0.5, color = "blue", size = 0.5) +
    geom_jitter(data = dados_confusao_teste_todos_falso_positivo , aes(x = valor, y = 0), alpha = 0.5, color = "red", size = 0.5) +
    facet_wrap(~atributo, scales = "free") +
    theme_minimal() +
    theme(
        aspect.ratio = .4,
        strip.text = element_text(size = 5),
        strip.background.y = element_rect(size = c(0.1,0.1)),
        strip.switch.pad.grid = unit(.01,"cm"),
        strip.switch.pad.wrap = unit(.01,"cm"),
        axis.text =  element_text(size = 4) ,
        axis.text.y = element_blank()
        
        )

COMUNICANDO OS RESULTADOS

Criando visualizações interativas com Shiny

Shiny é um pacote R que permite a criação de aplicações WEB interativas que ativam código em R, com as possibilidades de formatação de tabelas, geração de gráficos, execução de modelos etc que as bibliotecas do R disponibilizam.

Uma aplicação Shiny mínima

Os códigos desta seção que fala do Shiny não são executados diretamente, apenas são mostrados.

Uma aplicação Shiny mínima tem 4 itens:

Carregamento da biblioteca shiny
Definição da interface do usuário, ou seja, a página HTML com a qual o usuário vai interagir. Essa interface é definida com o uso da função fluidPage()
Definição do comportamento da aplicação com a função server()
Execução da função shinyApp(ui, server) para construir e iniciar uma aplicação.

library(shiny)
ui <- fluidPage(
  "Hello, world!"
)
server <- function(input, output, session) {
}
shinyApp(ui, server)

Uma aplicação simples funcional

Abaixo uma aplicação Shiny simples

Na função ui() do aplicativo shiny, definimos a interface.

No exemplo abaixo, além dos títulos, há um select e espaços para um gráfico e duas tabelas

Na função server() do aplicativo shiny, definimos o comportamento

No exemplo abaixo, definimos como o gráfico e as tabelas são gerados a partir do conteúdo do select contido na interface.

Repare que há duas tabelas, cada uma gerada com uma biblioteca diferente.

A biblioteca reactable permite que tenhamos mais controle sobre a estética

A biblioteca rhandontable parece mais com uma tabela de Excel, permitindo inclusive copiar e colar os números.

selectPaises <-     
    selectInput(
        "pais", 
        label = "País", 
        choices = gapminder$country,
        multiple = TRUE 
    )


ui <- fluidPage(
    h1("Dashboard Gapminder"),
    hr(),
    h3("Filtros:"),
    selectPaises,
    h3("Gráfico:"),
    plotOutput("grafico"),
    h3("Tabela bonita:"),
    reactableOutput("tabela_reactable"),
    h3("Tabela pros viciados:"),
    rHandsontableOutput("tabela_rhanson")
)

server <- function(input, output, session) {

    output$grafico <- renderPlot({
        gapminder %>% 
            filter(country %in% input$pais) %>% 
            ggplot() +
            geom_line(aes(x = year, y = gdpPercap, color = country )) +
            geom_point(aes(x = year, y = gdpPercap, color = country )) +
            theme_light()
        
    }) 

    
    output$tabela_reactable <- renderReactable({
        
        gapminder %>% 
            filter(country %in% input$pais) %>%
            select(
                country,
                year,
                gdpPercap
            ) %>% 
            pivot_wider(
                names_from = year,
                values_from = gdpPercap
            ) %>% 
            reactable(
                defaultColDef = colDef(
                    format = colFormat(digits = 0, separators = TRUE) 
                )
            )
            
    })
        
    
    output$tabela_rhanson <- renderRHandsontable({
        gapminder %>% 
            filter(country %in% input$pais) %>%
            select(
                country,
                year,
                gdpPercap
            ) %>% 
            pivot_wider(
                names_from = year,
                values_from = gdpPercap
            ) %>% 
            mutate_at(
                vars(matches("[0-9]{4}")),
                ~round(x = .x, digits = 0)
            ) %>% 
            rhandsontable(
                readOnly = TRUE
            ) %>% 
            hot_cols(
                format = "0,000",
                language = "pt-BR"
            )
        
    })
    
    
}


shinyApp(ui, server)

Programação reativa

Diferentemente da programação feita em scripts, a execução no Shiny não segue a ordem dos comandos. É o conceito da programação reativa

A execução das funções do tipo render (que já vimos), reactive, observe e observeEvent é preguiçosa. Só acontece quando os inputs são alterados.

No caso do nosso aplicativo simples, apenas quando o input “pais” é alterado, as funções render relativas aos outputs (gráfico e tabela) são executadas.

Produtores e consumidores

Temos 2 papéis para os objetos: produtores e consumidores

Os inputs são produtores e os outputs são consumidores.

Vamos começar a ver agora objetos que funcionam tanto como produtores como consumidores.

As reactive expressions do shiny são objetos que assumem os dois papéis.

Elas possibilitam que cálculos e tratamentos que são comuns a várias saídas sejam executados apenas uma vez.

Usando reactive para evitar duplicação

Na aplicação que desenhamos anteriormente, o select de país filtra da mesma forma os dados que são mostrados no gráfico e nas tabelas.

Podemos inserir um reactive que vai realizar esse tratamento que é comum a todas as saídas de uma só vez.

Assim evitamos mais facilmente, também, duplicar código.

Reactive expressions

Esta é a sintaxe das reactive expressions.

selectPaises <-     
    selectInput(
        "pais", 
        label = "País", 
        choices = gapminder$country,
        multiple = TRUE 
    )


ui <- fluidPage(
    h1("Dashboard Gapminder"),
    hr(),
    h3("Filtros:"),
    selectPaises,
    h3("Gráfico:"),
    plotOutput("grafico"),
    h3("Tabela bonita:"),
    reactableOutput("tabela_reactable"),
    h3("Tabela pros viciados:"),
    rHandsontableOutput("tabela_rhanson")
)

server <- function(input, output, session) {
    
    
    
    dados_filtrados <- reactive({
        
        gapminder %>% 
            filter(country %in% input$pais)
        
    })
    
    
    output$grafico <- renderPlot({
            dados_filtrados() %>% 
            ggplot() +
            geom_line(aes(x = year, y = gdpPercap, color = country )) +
            geom_point(aes(x = year, y = gdpPercap, color = country )) +
            theme_light()
        
    }) 
    
    
    output$tabela_reactable <- renderReactable({
        
        dados_filtrados() %>% 
            select(
                country,
                year,
                gdpPercap
            ) %>% 
            pivot_wider(
                names_from = year,
                values_from = gdpPercap
            ) %>% 
            reactable(
                defaultColDef = colDef(
                    format = colFormat(digits = 0, separators = TRUE) 
                )
            )
        
    })
    
    
    output$tabela_rhanson <- renderRHandsontable({

        dados_filtrados() %>% 
            select(
                country,
                year,
                gdpPercap
            ) %>% 
            pivot_wider(
                names_from = year,
                values_from = gdpPercap
            ) %>% 
            mutate_at(
                vars(matches("[0-9]{4}")),
                ~round(x = .x, digits = 0)
            ) %>% 
            rhandsontable(
                readOnly = TRUE
            ) %>% 
            hot_cols(
                format = "0,000",
                language = "pt-BR"
            )
        
    })
    
    
}


shinyApp(ui, server)

Dashboard exemplo

Este dashboard de exemplo tem algumas funcionalidades interessantes que vamos explorar nos próximos slides

Dashboard GD

library(shiny)
library(tidyverse)
library(reactable)
library(sf)
library(waiter)
library(ggmap)
library(RColorBrewer)
library(scales)

#### LEITURA DE DADOS ####

gd <- read_rds("dados/gd/gd.rds") %>% 
    mutate(
        unidade = 1
    )

#### VARIÁVEIS DE ESTADO GLOBAIS ####

ultimos_pontos <- gd
ultimo_grafico <- NA
ultimo_n_amostra <- 0

#### OPÇÕES DOS SELECTS ####

agrupadores <- c(
    "Setor econômico" = "USER_DscCl", 
    "Tarifa" = "USER_DscGr",  
    "Tipo de geração" = "USER_DscMo", 
    "Município" =  "USER_NomMu", 
    "Região" = "USER_NomRe", 
    "UF" = "USER_SigUF",
    "Tipo de fonte" = "USER_SigTi", 
    "Combustível" = "USER_DscCo",
    "Agente" = "USER_SigAg"
)


somas <- c(
    
    "Quantidade" = "unidade",
    "Potência (MW)" = "USER_MdaPo"
)

#### ELEMENTOS DE UI ####


mapa <- plotOutput(
    "mapa", 
    height = "600px",
    brush = brushOpts(
        id = "brush_mapa",
        #resetOnNew = TRUE,
        opacity = 0.2,
        stroke = "black",
        fill = "white",
        clip = FALSE,
        delay = 5000
    )
    
)


selectUFs <- selectInput(
    "ufs",
    "UFs:",
    choices = unique(gd$USER_SigUF) %>% sort(),
    multiple = TRUE
)

slider_n_amostra <- sliderInput(
    "n_amostra",
    label = "% amostrado no gráfico",
    min = 10,
    max = 100,
    step = 10,
    value = 10,
    post = "%"
    
)

slider_tamanho_ponto <- sliderInput(
    "tamanho_ponto",
    "Tamanho do ponto",
    min = 0,
    max = 2,
    step = 0.05,
    value = 0.05
    
)


agrupador <- selectInput(
    "campo_grupo",
    "Agrupar por:",
    choices = agrupadores
)


soma_por <- selectInput(
    "campo_soma",
    "Somar por:",
    choices = somas
)


tabela_count <- reactableOutput("info")

grafico_barras <- plotOutput("grafico_barras")

#### DIAGRAMAÇÃO DA UI ####


ui <- fluidPage(
    theme = "bootstrap.css",
    use_waiter(),
    titlePanel(
        fluidRow(
            column(width = 1, img(height = 40, src = "logo-epe-azul-15-anos.gif")),
            column(
                width = 11,
                tags$div(style = "height:10px"),
                h3("Dashboard Geração Distribuída")
                
            )
        )    
    ),
    tags$hr(),
    sidebarLayout(
        sidebarPanel(
            h4("Filtros" %>% tags$b()),
            wellPanel(
                selectUFs 
            ),
            h4("Configurações do gráfico" %>% tags$b()),
            wellPanel(
                agrupador,
                soma_por,
                slider_n_amostra,
                slider_tamanho_ponto
            ),
            width = 3
        ),
        mainPanel(
            width = 9,
            splitLayout(
                cellWidths =  c("60%","40%"),
                div(style = 'overflow: hidden',mapa),
                div(style = 'overflow: hidden',
                    verticalLayout(
                        grafico_barras,
                        tabela_count
                    )
                )
            )
        ) 
    )
)




server <- function(input, output, session) {
    
    #### MAPA ####
    
    
    output$mapa <- renderPlot({
        
        
        pontos_selecionados <- pontos()

        if(isTRUE(all_equal(ultimos_pontos, pontos_selecionados)) & ultimo_n_amostra == input$n_amostra)
        {
            resposta <- ultimo_grafico
        }
        else
        {
            print("atualiza")

            xmin <- min(pontos_selecionados$X)
            ymin <- min(pontos_selecionados$Y)
            xmax <- max(pontos_selecionados$X)
            ymax <- max(pontos_selecionados$Y)
            
            margem_x <- min((ymax - ymin)/7,2)
            margem_y <- min((xmax - xmin)/7,2)
            
            xmin <- xmin - margem_x 
            ymin <- ymin - margem_y
            xmax <- xmax + margem_x 
            ymax <- ymax + margem_y
            
            
    
            if (is.na(xmin)){
                local <- c(-67.85551, -31.76278, -34.80921, -1.198088)
            }
            else{
                local <- c(xmin, ymin, xmax, ymax)    
            }
            
            map = get_map(local,source = "stamen", maptype = "toner-lite") 
            
            
            resposta <- ggmap(ggmap = map) +
                stat_density_2d(
                    aes(x = X, y = Y, fill = ..level.. ), 
                    bins = 30,
                    geom = "polygon", 
                    data = pontos_selecionados ,
                    alpha = .1
                ) +
                geom_point(
                    data = pontos_selecionados, 
                    aes(
                        x = X,
                        y = Y
                    ),
                    alpha = 0.01,
                    color = "darkblue",
                    size = input$tamanho_ponto
                ) +
                scale_fill_gradientn(colors = brewer.pal(7, "YlOrRd")) +  
                theme_inset()
            
            ultimo_grafico <<- resposta 
            ultimos_pontos <<- pontos_selecionados
            ultimo_n_amostra <<- input$n_amostra 
        
        }

        resposta 
        
    })

    #### PONTOS DO MAPA SELECIONADOS ####

        
    pontos <- reactive({
        
        
        waiter <- Waiter$new(id = "mapa")$show()

        resposta <- brushedPoints(gd_sample(), input$brush_mapa, xvar = "X", yvar = "Y") 
        

        if (nrow(resposta) == 0){
            resposta <- gd_sample()
        }
        

        resposta
    })

    #### DADOS SELECIONADOS PELO MAPA####
    
        
    dados <- reactive({

        resposta <- brushedPoints(gd_filtro(), input$brush_mapa, xvar = "X", yvar = "Y") 
        
        if (nrow(resposta) == 0){
            resposta <- gd_filtro()
        }

        campo_grupo <- names(agrupadores)[agrupadores == input$campo_grupo]
        campo_soma <- names(somas)[somas == input$campo_soma]
        
        resposta %>% 
            group_by_at(input$campo_grupo) %>% 
            summarise_at(input$campo_soma, sum) %>% 
            ungroup() %>% 
            mutate(
                Parcela = .data[[input$campo_soma]]/sum(.data[[input$campo_soma]])
            ) %>% 
            arrange(desc(.data[[input$campo_soma]])) %>% 
            rename(
                !!campo_grupo := 1,
                !!campo_soma := 2
            )
        
    })

    #### TABELA ####
    
    
    output$info <- renderReactable({
        

        
        dados() %>% 
            reactable(
                pagination = FALSE,
                defaultColDef = colDef(
                    format = colFormat(
                        digits = 0,
                        separators = TRUE
                    )
                ),
                columns = list(
                    Parcela = colDef(
                        format = colFormat(
                            percent = TRUE,
                            digits = 0
                        )
                    )
                )
                    
            )
        
    })

    #### GRÁFICO DE BARRAS HORIZONTAIS ####
    
    
    output$grafico_barras <- renderPlot({

        campo_grupo <- names(agrupadores)[agrupadores == input$campo_grupo]
        campo_soma <- names(somas)[somas == input$campo_soma]

        dados() %>% 
            mutate(
                !!campo_grupo := 
                    fct_lump(
                        f = .data[[campo_grupo]], 
                        prop = 0.05, 
                        w = .data[[campo_soma]] ,
                        other_level = "Outros"
                    )
            ) %>% 
            mutate(
                !!campo_grupo :=
                    fct_reorder(
                        .f = .data[[campo_grupo]],
                        .x = .data[[campo_soma]],
                        .fun = sum
                    )
            ) %>% 
            ggplot() +
            geom_col(
                aes(
                    x = .data[[campo_grupo]],
                    y = .data[[campo_soma]],
                    fill = .data[[campo_grupo]]
                )
            ) +
            labs(
                y = campo_soma,
                x = campo_grupo,
                fill = campo_grupo
            ) +
            coord_flip() +
            theme_minimal() +
            scale_y_continuous(
                labels = number_format(
                    big.mark = ".",
                    accuracy = 1
                )
            ) +
            theme(
                legend.position =  "top"
            ) +
            NULL
        
        
        
    })
    
    

    gd_filtro <- reactive({
        
        gd %>%
            filter(USER_SigUF %in% input$ufs | length(input$ufs) == 0 )
        
    })
    
    gd_sample <- reactive({
        gd_filtro() %>% 
            sample_frac(input$n_amostra/100)
    })


}


shinyApp(ui, server)

Mapa interativo

O usuário pode selecionar áreas do mapa.

A seleção dessa área leva a uma seleção de pontos

Essa opção é habilitada com brushOpts e os pontos selecionados são acessados com brushedPoints

Este tipo de interação pode ser usado com qualquer tipo de gráfico

Spinner de espera

Esta aplicação lida com mais de 100.000 dados a respeito de instalações de geração distribuída.

Portanto, o gráfico demora um pouco para ser plotado.

Para que o usuário saiba que a aplicação está calculando alguma coisa e não travada, usamos um waiter

Uso de váriáveis globais de estado

Algumas variáveis de estado são criadas para ajudar a controlar a aplicação.

Elas são atribuídas com o operador <<-

Programando com Tidy evaluation

O Tidy evaluation deixa que usemos os nomes dos campos sem aspas.

Em aplicações Shiny exploratórias, é comum que os campos que participam dos agrupamentos ou dos cálculos sejam escolhidos pelo usuário. Portanto, estes campos estarão no comteúdo de variáveis.

Para que as funções do tidyverse usem o conteúdo da variável como o campo a ser considerado (e não o próprio nome da variável), temos que lançar mão das funções group_by_at, summarise_at, .data[[variável]] e o operador !!.

Diagramação da página

A diagramação da página é toda feita com as funções de diagramação do Shiny.

Entretanto, cada objeto é definido anteriormente, o que é diferente do template que é gerado pelo RStudio quando um template Shiny é gerado.

Isso parece besteira, mas separar a criação dos objetos da diagramação destes objetos facilita muito uma reorganizaçào da diagramação e causa menos problemas relacionados a erros de sintaxe como falta de vírgulas e parênteses

Criando relatórios com R Markdown

Livros

Referência Bookdown

tidyverts

Os pacotes integrantes do conjunto tidyverts tratam as séries temporais do jeito tidy

-Funções básicas para o tsibble (tibble temporal) tsibble

-Funções pra extrações de features e estatísticas feasts

-Funções para projeção fable

tsibbles

***tsibble%%% é uma classe derivada da tibbl que tem funcionalidades para séries temporais

ipca <- BETS::BETSget(433)

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

A função ***as_tsibble()%%% cria um tibble a partir de vários tipos de objetos que armazenam séries temporais

O tsibble representa bem séries temporais regulares e reconhece o intervalo entre os valores

tsbl_ipca <- as_tsibble(ipca) %>% 
  mutate(value = value/100)


print(tsbl_ipca)

## # A tsibble: 481 x 2 [1M]
##       index  value
##       <mth>  <dbl>
##  1 1980 jan 0.0662
##  2 1980 fev 0.0462
##  3 1980 mar 0.0604
##  4 1980 abr 0.0529
##  5 1980 mai 0.057 
##  6 1980 jun 0.0531
##  7 1980 jul 0.0555
##  8 1980 ago 0.0495
##  9 1980 set 0.0423
## 10 1980 out 0.0948
## # ... with 471 more rows

tsibbles com várias séries a partir de um tibble

No código abaixo calculamos o número de instalações de Geração Distribuída por UF as partir dos dados a respeito das instalações.

Repare que a função ***yearmonth()### cria uma data mensal que é específica para períodos de um mês

dados_gd <- read_rds("dados/gd/gd.rds")

meses <- tibble(mes = seq.Date(from = make_date(2009,1,1), to = make_date(2020,3,1), by = "month" ) %>%  yearmonth())


dados_gd_cum <- dados_gd %>% 
    rename(
        uf = USER_SigUF
    ) %>% 
    crossing(meses) %>% 
    filter(DthConexao < mes) %>%
    group_by(
        uf,
        mes
    ) %>% 
    summarise(
        n = n(),
        potencia = sum(USER_MdaPo)
    )

A função as_tibble aceita um campo como referência da data, via parâmetro index, e um campo como chave de cada série contida no tibble: key

dados_gd_tsbl <- as_tsibble(dados_gd_cum, key = uf, index = mes ) 


dados_gd_tsbl

## # A tsibble: 2.179 x 4 [1M]
## # Key:       uf [28]
## # Groups:    uf [28]
##    uf         mes     n potencia
##    <chr>    <mth> <int>    <dbl>
##  1 AC    2015 mai     1      4  
##  2 AC    2015 jun     1      4  
##  3 AC    2015 jul     2      6  
##  4 AC    2015 ago     3     38.5
##  5 AC    2015 set     3     38.5
##  6 AC    2015 out     3     38.5
##  7 AC    2015 nov     3     38.5
##  8 AC    2015 dez     3     38.5
##  9 AC    2016 jan     3     38.5
## 10 AC    2016 fev     3     38.5
## # ... with 2.169 more rows

Mudando a periodicidade

Podemos mudar a periodicidade da série indexando por uma transformação do índice atual, usando group_by_key%%% e index_by%%%.

No caso, queremos as informações no final do ano, portanto vamos usar ***last()%%%.

gd_anual <- dados_gd_tsbl %>% 
    group_by_key() %>% 
    index_by(ano = year(mes) ) %>% 
    summarise(
        n = last(n),
        potencia = last(potencia)
    )

gd_trimestre <- dados_gd_tsbl %>% 
    group_by_key() %>% 
    index_by(trimestre = yearquarter(mes)) %>% 
    summarise(
        n = last(n), 
        potencia = last(potencia)
    )


gd_anual

## # A tsibble: 214 x 4 [1Y]
## # Key:       uf [28]
##    uf      ano     n potencia
##    <chr> <dbl> <int>    <dbl>
##  1 AC     2015     3     38.5
##  2 AC     2016     4     40.5
##  3 AC     2017    12    173. 
##  4 AC     2018    38    630. 
##  5 AC     2019   134   1968. 
##  6 AC     2020   172   2353. 
##  7 AL     2015     3     57.1
##  8 AL     2016    18    270. 
##  9 AL     2017    77   1195. 
## 10 AL     2018   255   3203. 
## # ... with 204 more rows

gd_trimestre

## # A tsibble: 737 x 4 [1Q]
## # Key:       uf [28]
##    uf    trimestre     n potencia
##    <chr>     <qtr> <int>    <dbl>
##  1 AC      2015 Q2     1      4  
##  2 AC      2015 Q3     3     38.5
##  3 AC      2015 Q4     3     38.5
##  4 AC      2016 Q1     3     38.5
##  5 AC      2016 Q2     3     38.5
##  6 AC      2016 Q3     4     40.5
##  7 AC      2016 Q4     4     40.5
##  8 AC      2017 Q1     4     40.5
##  9 AC      2017 Q2     7     75.8
## 10 AC      2017 Q3     9    121. 
## # ... with 727 more rows

Plotando uma série

A função ***autoplot()### cuida das principais tarefas para plotar uma série

tsbl_ipca %>% 
  filter(index > make_date(1994,7,1)) %>% 
  autoplot(value, color = "darkblue", size = 1) +
  scale_y_continuous(labels = percent_format(accuracy = 1) ) + 
  scale_x_date(date_breaks =  "1 year", labels = year) +
  labs(
    x = "",
    y = "IPCA"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90)
  )

Plot sazonal

A função gg_season ajuda a fazer um gráfico mostrando os períodos sazonais para comparação

tsbl_ipca %>% 
  filter(year(index) > 2009) %>% 
  gg_season(value, period = "year", labels = "both") +
  scale_y_continuous(
    labels = percent_format(accuracy = 0.1), 
    breaks = seq(0.000, by = 0.001,to =  0.015) ) + 
  scale_x_date(
    date_labels = "%b",
    date_breaks = "1 month"
  ) +
  labs(
    x = "",
    y = "IPCA"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90)
  )

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

tsbl_ipca %>% 
  filter(year(index) > 1995) %>% 
  gg_season(value, period = "year") +
  scale_y_continuous(
    labels = percent_format(accuracy = 0.1), 
    breaks = seq(0.000, by = 0.002,to =  0.03) ) + 
  scale_x_date(
    date_labels = "%b",
    date_breaks = "1 month"
  ) +
  labs(
    x = "",
    y = "IPCA"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90)
  )

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

buscas <- gtrends("xvideos", geo = "BR", time = "today 1-m", onlyInterest = TRUE)

buscas$interest_over_time

##          date hits geo      time keyword gprop category
## 1  2020-02-04   86  BR today 1-m xvideos   web        0
## 2  2020-02-05   89  BR today 1-m xvideos   web        0
## 3  2020-02-06   87  BR today 1-m xvideos   web        0
## 4  2020-02-07   90  BR today 1-m xvideos   web        0
## 5  2020-02-08  100  BR today 1-m xvideos   web        0
## 6  2020-02-09   98  BR today 1-m xvideos   web        0
## 7  2020-02-10   87  BR today 1-m xvideos   web        0
## 8  2020-02-11   86  BR today 1-m xvideos   web        0
## 9  2020-02-12   83  BR today 1-m xvideos   web        0
## 10 2020-02-13   84  BR today 1-m xvideos   web        0
## 11 2020-02-14   86  BR today 1-m xvideos   web        0
## 12 2020-02-15   98  BR today 1-m xvideos   web        0
## 13 2020-02-16   94  BR today 1-m xvideos   web        0
## 14 2020-02-17   85  BR today 1-m xvideos   web        0
## 15 2020-02-18   84  BR today 1-m xvideos   web        0
## 16 2020-02-19   83  BR today 1-m xvideos   web        0
## 17 2020-02-20   83  BR today 1-m xvideos   web        0
## 18 2020-02-21   89  BR today 1-m xvideos   web        0
## 19 2020-02-22   96  BR today 1-m xvideos   web        0
## 20 2020-02-23   98  BR today 1-m xvideos   web        0
## 21 2020-02-24   97  BR today 1-m xvideos   web        0
## 22 2020-02-25   98  BR today 1-m xvideos   web        0
## 23 2020-02-26   87  BR today 1-m xvideos   web        0
## 24 2020-02-27   81  BR today 1-m xvideos   web        0
## 25 2020-02-28   85  BR today 1-m xvideos   web        0
## 26 2020-02-29   92  BR today 1-m xvideos   web        0
## 27 2020-03-01   93  BR today 1-m xvideos   web        0
## 28 2020-03-02   83  BR today 1-m xvideos   web        0

DESENVOLVIMENTO DE PACOTES

Filosofia

Tudo que pode ser automatizado deve ser automatizado devtools

Preparando o ambiente

Pacotes a instalar para começar:

install.packages(c("devtools", "roxygen2", "testthat", "knitr"))

Para checar se o ambiente está preparado para criar pacotes:

devtools::has_devel()

Criando a estrutura

New Project -> New Directory -> R Package

The main difference between library() and require() is what happens if a package isn’t found. While library() throws an error, require() prints a warning message and returns FALSE. In practice, this distinction isn’t important because when building a package you should NEVER use either inside a package. See package dependencies for what you should do instead.